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■ Abstract 

,—1 ■ In this paper, we study the evolution of a relativistic, superhorizon-sized 

'i> ', void embedded in a Friedmann-Robertson- Walker universe. We numerically 

' solve the spherically symmetric general relativistic equations in comoving, 

^ I synchronous coordinates. Initially, the fluid inside the void is taken to be 

^ ' homogeneous and nonexpanding. In a radiation-dominated universe, we find 

Q I that radiation diffuses into the void at approximately the speed of light as a 

' strong shock — the void collapses. We also find the surprising result that the 
cosmic collapse time (the I'^^-crossing time) is much smaller than previously 

^ ■ thought, because it depends not only on the radius of the void, but also on 

P-i! the ratio of the temperature inside the void to that outside. If the ratio 

p • of the initial void radius to the outside Hubble radius is less than the ratio 

Qj [ of the outside temperature to that inside, then the collapse occurs in less 

• than the outside Hubble time. Thus, superhorizon-sized relativistic voids 

^ I may thermalize and homogenize relatively quickly. These new simulations 

• '-j ■ revise the current picture of superhorizon-sized void evolution after first-order 

r> ' inflation. 
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I. Introduction 

The Big Bang model predicts tlie evolution of a liomogeneous and isotropic universe. 
The confirmation of its predictions (e.g. the ?)K Planck spectrum of the microwave back- 
ground, the primordial abundances of ^ife, ^He, D, and "^Li from Big Bang Nucleosyn- 
thesis, and the number of light neutrino species) is stunning. One is led to question why 
the universe can be described so well by a homogeneous and isotropic model. 

Inflation can provide the answer to this question. It occurs when a scalar field, a has 
non-zero potential energy V{a) which dominates the energy density of the universe.EJ'EI The 
key ingredients to all models are that the universe expands superluminally during inflation 
and that there is massive entropy generation afterward. If the universe increases at least 
10^^ times its original size, the flatness and horizon problems are solved; the universe is 
dynamically driven to homogeneity and isotropy. The most interesting class of models is 
first-order inflation. Here, the scalar field is trapped in the false vacuum state of a strongly 
first-order potential. Bubbles of true vacuum are nucleated at different spacetime points, 
and the end of inflation occurs when the universe is fillejd with true-vacuum bubbles of 
varying sizes. The original model, Guth's "old inflation" p does not work because it fails 
to percolate (fill) the universe with true-vacuum bubbles. More recent models of first- 
rder inflation which modify the gravitational or particle sector (e.g. "extended inflation" 
) are promising as early-universe inflationary scenarios. Because the universe expands 
as a power-law in time rather than exponentially, percolation is guaranteed to occur. 

The end of inflation occurs when true-vacuum bubbles of different sizes fill all of 
space. A confusing mess of scalar field dynamics then occurs as the bubbles collide.cl 
Because all the energy is contained in the bubble walls, reheat occurs when the cr-field 
gradient energy is converted into locally thermal radiation. After reheat, the standard 
homogeneous and isotropic Big Bang model describes the evolution of that part of the 
universe that had contained horizon-sized bubbles, since these would have been ther- 
malized during reheat. The very large superhorizon-sized bubbles (which were nucleated 
early on during inflation), however, have traditionally been a problem.B After the a-field 
in the bubble wall decays to relativistic particles, a nearly empty void is formed. Since 
these voids are much larger than the Hubble radius outside the void, it has been thought 
that the inside of the largest superhorizon-sized voids would remain empty until after 
recombination, thus producmg unobserved temperature fluctuations of order unity in 
the microwave background.&LrQ The longevity of these voids follaws from assuming that 
a superhorizon-sized void expands conformally with spacetime. &□ The earliest time at 
which thermalization could occur would then be when the Hubble radius outside the 
void is of order the size of the void, since this is the expected f^^-crossing time (i.e. the 
time for photons originally in the void wall to reach the origin). If the void has comoving 
size ro, this occurs in time At = t — U ^ /7^t(ti)(c^^-R„aii(ti)/-ff,^l(ti))^, where the initial 
void radius is -Rwaii(''^i) = roa(ti), cH^^^^iti) is the Hubble radius outside the void, a(t) is 
the cosmic scale factor and t[ is the cosmic time after reheating. Other authors suggest 
that the void would fill in with radiation, but this is estimated to occur on similar time 
scales.B Because of this "big bubble problem" , first-order inffationiaHodels are fine-tuned 
in order to keep the production of large bubbles at a minimum. oBa 

Motivated by the first-order (e.g. extended) inflation "big-bubble problem", we 
present numerical studies of the evolution of a superhorizon-sized general relativistic void 
embedded in a Friedmann-Robertson- Walker universe. Altliqugh pressureless and thin- 
shell superhorizon-sized voids have been studied in the past,liy'li3 this is the first study of 
general relativistic voids with pressure and of arbitrary size and void wall structure. We 
emphasize that our results are not dependent on any particular inflation model. These 
new simulations show that opposite sides of a superhorizon-sized relativistic void interact 
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in a very short cosmic time, thereby suggesting that these voids can also thermahze and 
homogenize on short time scales. 

In Section II, we discuss the general relativistic spherically symmetric metric in La- 
grangian gauge and synchronous coordinates, and present the equations to be solved 
numerically. In Section III, the initial conditions, boundary conditions and numerical 
techniques used are described. In addition, remarks about the deceleration of a void wall 
are made. In Section IV, we derive the surprising result that the I'^^-crossing time of 
a superhorizon-sized general relativistic void can be vanishingly short. Section V con- 
tains numerical solutions for pressureless dust and comparisons to exact Tolman-Bondi 
solutions. A test of the code for the relativistic Friedmann-Robertson- Walker solution is 
presented in Section VI. The numerical evolution of nonrelativistic, special relativistic, 
and general relativistic voids with pressure is examined in Section VII. Here it is found 
that a relativistic superhorizon-sized void collapses in the form of a strong shock moving 
at the speed of light. Thus, the collapse time is approximately the l^*-crossing time. 
Finally, Section VIII contains a discussion of the results. 

II. Spherically Symmetric General Relativistic Fluids 
A. The (Lagrangian) Metric 

The most general spherically symmetric metric isli^ 

ds"^ = cH{t, r)dt^ + a{t, r)drdt + h{t, r)dr'^ + k{t, r)d^? , (2.1) 

where dVt^ = dO"^ + sin^ Odtp"^ and c is the speed of light. To this metric we can apply 
general transformations of the type t = fi{t',r') and r = /2(t',r') without altering the 
spherical symm etry. If we perform the necessary transformations to eliminate the drdt 
term, then Eq. can be written as 

ds^ = -c2$2(t, r)dt'' + A\t, r)dr^ + R^{t, r)dn^, (2.2) 

where we have dropped the primes, and where we require our coordinates to be comoving 
with the fluid. At time t, the metric function R{t, r) is the Eulerian distance that a fluid 
shell labeled by r is located from the center of coordinates. More precisely, 27ri?(t, r) is 
the spacelike circumference of a sphere centered on the origin which contains all particles 
with comoving coordinate r. We have thus chosen the Lagrangian gauge with synchronous 
coordinates (Gaussian normal coordinates). Transformations of the form t = f{t) and 
f = g{r) can still be made. This metric was first used to study the eeneral relativistic 
collapse of stars to black holes or neutron stars during supernovae.EJa 

It is very important in numerical general relativity to carefully choose the appropri- 
ate gauge and coordinates to best match the physical problem to be studied. We have 
thus specifically chosen the Lagrangian gauge and synchronous coordinates to evolve a 
general relativistic void embedded in an expanding Friedmann-Robertson- Walker (FRW) 
universe. In the Lagrangian gauge, we gain maximal coverage of the fluid in a numerical 
scheme. This is important since the void is embedded in an expanding universe, so that 
we would continuously lose mass shells in a Eulerian scheme. In addition, the final re- 
sults are much more easily relatable to our own approximately FRW homogeneous and 
isotropic universe. We note that asynchronous coordinates could instead be used with 
the Lagrangian gauge. However, here time and space are mixed up so that the region 
outside a void is no longer spatially homogeneous — the necessary initial conditions are 
not obvious and would need to be determined using comoving synchronous coordinates 
to initially set up and evolve the void. Asynchronous slicing of space-time (e.g. polar 
slicingt^) is usually used numerically to study the collapse of a star to a black hole in 
a flat non-expanding universe. These coordinates are necessary to study the physically 
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interesting mass zones outside the apparent horizon after a mass shell has crossed this 
horizon. (In synchronous coordinates, once a mass shell crosses the apparent horizon, 
numerical integration stops so that the evolution of the mass shells outside this horizon 
cannot be studied). Since we are evolving an underdense region here and do not expect 
apparent horizons to form, we do not need to resort to these coordinates. 

Because we wish to embed a void in a FRW universe, we now relate the metric 
functions from Eqn(|2.2|) to those from the familiar FRW metric: 



ds^ = -cV + a{tf (^^^ + r-dn-) , (2.3) 

where a{t) is the cosmic scale factor, and is —1, or 1 for negative, zero and positive 
spatial curvature, respectively. The solutions describe a universe which is homogeneous 
and isotropic on each time slice.i For the FRW metric then, $ = 1, A = a{t) / vT^^XcFV^ 
and the Eulerian distance is R{t) = ra{t). 

In this paper, the particles are assumed to be everywhere in local thermal equilibrium, 
so that we can describe them as a fluid. The stress-energy tensor for a viscous fluid with 
energy density p and pressure p (measured in the frame of the fluid) is 

= pu'^u^ + pP''^ - 2r](T'^^ - 3CeP°^ (2.4) 

where is the fluid 4— velocity, P"^ = + g""^ is the projection operator, a"^ = 
1/2 (VpM"P'''^+V;,M^P^")-l/3 eP"^ is the shear viscosity tensor, 6 = V^m^ is the fluid 
expansion coefficient, an d and ( are arbitrary functions of r and t. If r/ = (, = 0, the fluid 
is non-viscous. For Eqn (|2.2|) with comoving fluid 4-velocity = (— 0, 0, 0), a/ = (5 

and a/ = a^,^ = -1/2 /?, where {3 = 2/3 (A/A - U/R), and 6 = A/A + 2U/R. (For 
the FRW metric, (3 = and G = 3a/ a). The stress tensor is diagonal, with components 

Tt' = -p, T/ = {p-Ce)-2r](3, and T/ = = (p-CQ) +r//5. We will employ a scalar 
artificial viscosity, Q, so that rj = and (Q = —Q, where Q will be significantly non-zero 
only in areas of steep "velocity" gr-adients, as in a shock. This is essential for stabilizing 
numerical shocks, as is well known.Ej This viscosity will dissipate enough energy on small 
scales so that the numerical solution approaches the exact solution in the limit that the 
grid spacing approaches zero. 

B. Fluids composed of massless particles 
If a fluid consists of (effectively) massless, locally thermalized particles, we can relate p 
to p through the equation of state p = pLu)- Fo'^ photons (or particles with mass p having 
local temperatures T :^ p), p = p/3. lO Setting TZj - gjTl = SttGnTJ, where TZj 
is the Ricci tensor and is Newton's constant, and using the conservation equations 
y^ji/io _ independent equations are found: Gq^ = To°, Gi^ = Ti^, Gq^ = 

To\ TO'^;^ = 0, and T^^;^ = 0. We define 

r = P'/A, (2.5) 
M' = Anc-^pR^R' (2.6) 

and U = ^~^{dR/dt), where the prime denotes differentiation with respect to r. The 
general relativistic equations can be written 

P = $f/ (2.7) 



We set c = 1 for the rest of this subsection 
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. ^[Gr^M , 4nG^{p + Q)R \ c'T^^p + Q)' 

- )~{p + p + Q)R' ^^-^^ 

M = -4tt{p + Q)R'^^U/c^ (2.9) 

p = -^(p + p + Q)^-^ (2-10) 

= -$J^±^ (2.11) 

P+P+Q ^ ' 

= l + (t//c)^-2G^M/(i?c2), (2.12) 

where the dot denotes differentiation with respect to t, and where we have included an 
alternate definition for F. There is also the auxiliary equation: A = A$f/'/i?. 

The quantity IJ = describes a particle's "velocity" as measured in the frame of 
the fluid, since for dr = d6 = dtp = 0, the infinitesimal proper time that each observer 
measures is dr = \/—ds'^ = c^dt. For Gat = 0, if a particle has velocity v, then 
F = 1/^1 — {v/cY and U = Tv (see Eqn ( p.9|) ); F and U represent the two non-trivial 

components of the 4- velocity of the fluid. If F ^ 1, then f//c ~ F and the fluid is moving 
at relativistic velocities relative to a stationary observer. 

We can rewrite Eqn ( p.6| ) in terms of the proper volume element d^V = AnR'^Adr = 
AT^R^R'dr/T. The "mass-energy" function M then becomes 

M{t,r)=c-'^ j d^VVp. (2.13) 

For Gtv = 0, because the volume along the radial direction is Lorentz-contracted by the 
factor 1 / \Jl — (f/c)^ = F and p is the energy density measured in the frame of the 
fluid, Fp is just the energy density of a fluid parcel as measured by a stationary observer. 
Therefore, M(r, t) is just the total "mass-energy" contained within comoving coordinate 
r at time t. 

To better understand these equations, we relate them to the FRW equations. Recall 
the FRW results following Eqn (|2.3| ): $ = 1, i? = ra and A = a/ \Jl — kr"^ /c^. From 



Eqn (|2.7|) , the "velocity" is U = R = rd. Since R' = a(t), F = — from 

Eqn (|2.5|) . For a spatially flat [k = 0) FRW model then, F = 1 even though the fluid at 
R is moving away from the origin with velocity U. Thus in general relativity, F in not 
the relativistic gamma-factor of the fluid with respect to the origin. Using the fact that 
M = Attc~'^pR^ /3, we see that Eqn ( I2.12| ) is just Friedmann's equation, = (a/a)^ 
SttGnc-^p/S - kc^a^ with H = $"V7a = U/R. 



In a spatially flat FRW universe then, U = c'^R^StiGnp/^ = ^J2GnM/R. If we 
set U = f/cRAv + UpEC: where ?7grav = IG^M/R is the gravitational velocity and 
UpEc is the peculiar velocity, F can be expressed in terms of these velocities: F^ = 
1 + {UpEc/cf + 2UpecUgkav/c\ 

For a non- viscous fluid with p = (7 — l)p, Eqn ( p. 11] ) (the conservation of momen- 
tum equation, T^i.^ = 0) can be integrated exactly to give ^(t,r2) = $(t, ri)[p(t, ri)/ 

p(t,r2)]^'''~^^^"' . We are interested in evolving a void embedded in a FRW homogeneous 
and isotropic universe, so that p' = p' = outside the void. We define $out and pout 
to be the spatially constant values of $ and p outside the void. (In what follows, the 
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subscripts "in" and "out" represent the spatially constant values inside and outside the 
evolving void region, respectively). Taking 7 = 4/3, we find 

$=$out(^y^'. (2.14) 

The fact that this equation can be solved exactly is important for calculating the 1*^*- 
crossing time for general relativistic voids, as will be discussed in Section IV. 

To gain some physical insight into $, we calculate the potential energy of a fiuid 
distribution. A comoving observer has dr = d6 = dijj = 0, so that the proper acceleration 

measured by this observer is = $^^f/. Since the force is radial, we can write ar = 
—d(l)/dR, where (p is the potential. Using Eqn( p.8|) and integrating, we can write the 
potential as the following sum: R) = 0grav(^, R) + 0fluid(^, -R), where 



5'GRAV = 4:7iGn 







R r -I 

M/R^ + c-^pR dR (2.15) 



^FLUID 



/ c^V^dp/{p + p) (2.16) 

JO 

for (5 = 0. For a void with p = p/3, we obtain the usual gravitational potential 
0GRAV — 47rG7vPin-R^/3 = Gj\fM/R inside the void, and 0grav — 47rG7vPout-R^/3 — 
G]\fM/R outside. If in addition we choose T{ti,R) = 1 initially, then 0FLUiD(^i) = 
ln[ p{ti, R)/ piniU) ]^/^ = -ln[ $(ti,i?)/$in(ti) ]. (Note that the contribution to the fiuid 
potential is zero inside, but (potentially much) greater than zero outside the void. Thus 
it can substantially increase the already large potential outside the void). Therefore, $ 
is proportional to the exponential of the fiuid potential, 0fluid- 

C. Fluids composed of massive particles 
Suppose instead we consider a fiuid which consists of particles of mass p and with 
arbitrary temperature T. Then, the total energy density of the fiuid is the mass energy 
density plus the internal energy density. Denoting the (proper) mass density {p times 
the number density) by n{t,r) and the internal energy per unit mass by e{t,r), 

p = c^n{l + e/c^). (2.17) 

Here we have traded one unknown for two because in general the pressure depends not 
only on p but also on n. For a fiuid composed of relativistic (nonrelativistic) particles, 
e/c^ > 1 (e/c^ < 1). If the fiuid obeys the perfect gas law, then its pressure is 

p={-f-l)ne. (2.18) 

For a highly relativistic species with 7 = 4/3, the energy density is three times the 
pressure: p = ne = p/ {'y — 1) = 3p, whereas for a highly nonrelativistic fiuid, the energy 
density is much larger than the pressure: p ~ nc^ = p/[(7 — l)e/c^] ^ p. We assurne 
that the total number of particles per comoving volume is constant: V ^{nu^) = 0. 
Using Eqn(P3|), this can be integrated to give 

/(r) = AimR^R'/T (2.19) 
= (2.20) 

where /(r) is an arbitrary function depending only on the coordinate r. Specifying /(r) 
completely fixes the arbitrariness of the metric functions under transformations in r (as 
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discussed after Eqn ( p.2[ )). This particular definition for / is necessary in order to write 
the difference schemes in a geometrical way that allows shocks and explosions to be 

numerically stable at the origin.li3 r~\r~\ 

We can now rewrite the full set of general relativistic equations ( (|2.7| )- (|2.12| )) agi^li^l 

R = (2.21) 

U = ^ ( ^^^ I + Q^^ ] 47rr$i?^(p + Q)' 

\ R'^ J wr"^ 

M = -A'n{p + Q)R^^U/c^ (2.23) 
. = _ 4xt(p + 9)(fl^t/)- p.26) 



nwc 



2 



(2.26) 



where F is given by E qn (|2.12|) and w = 1 + {e+p/n)/(? is the relativistic enthalpy. 



Equations ( 2.21|) -( p.26|) (along with the definitions for F and w given in the previous 
sentence) are the set used in the numerical codej^ 

When the kinetic energy of each particle is much less than its mass energy e/c^ <^ 1 
(or T/yU ^ 1), we obtain the nonrelativistic Lagrangian fluid equations.EiJ (They can also 
be obtained by setting — oo in Eqns ( p.21D -( p.26| )). For future reference, in this limit 
$ — > 1, F ^ 1, and w — > 1 so that U = Ris the fluid velocity, M{t, r) = r^/3 is the total 
mass within r and n = r^/ {R?R') is the mass density. 

The artificial viscosity used here is given by Equation ( |C21| ), and is generalized from 
the expression used by previous workers:B3o 

Q = k\{l + e/{rc^)){U'fdr^/r for U' < 

g = otherwise. (2.27) 

Consider the behavior of n and p in the limit that the entropy, S{t,r), within r is 
conserved. Using the thermodynamic relation TdS = d{p/n) +pd{l/n), we find e = 
pri/n^. Then, the pressure and internal energy for a shell labeled by r are related to their 
initial values hy p (x and e oc In the ultrarelativistic limit with 7 = 4/3, 

n oc p^/^ and e oc p^^^. In addition, if the fluid is a ploytrope (L£. isentropic), 5" = 0. 
Using TdS = d{p/n) +pd{l/n) again, we find the familiar resultE3 

(7-l)/7 

I • (2.28) 



e{t,ri) fn{t,ri)y fp{t,ri) 



e{t,r2) \n{t,r2)J \p{t,r2) J 

For ultrarelativistic fluids, p = ne and therefore e oc p(T~^)/'>' and n oc p^^"'. If 7 = 4/3, 

e oc p^/^ and n oc p^/^. This is the property of a relativistic fluid, since then p (xT^ and 
n oc T^. 



''For technical difficulties in the general relativistic case, we determine M via the M equation rather 
than the M' equation. In addition, we determine n via the n equation rather than through the analytical 
solution n — Tr^ /{R^R'), because too much intrinsic viscosity is introduced for special relativistic voids 
otherwise. 
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D. Fluid Deceleration in the Void Wall 

In this section we show that if the outward peculiar velocity of the wall of a general 
relativistic void is very large, then the deceleration of this wall can be enormous. This 
"damping force" is responsible for slowing down and collapsing a superhorizon-sized 
void. Without this fluid force, the void would expand (not collapse) from gravitational 
forces(see V). 

In section 11-B, we saw that for special relativistic fluids {Gn = 0), U is the net 
outward momentum per particle mass fi. Similarly, if in the general relativistic case 
^ ^GRAv — '^GnM/ R, the gravitational attraction inward is negligible and the 
dynamics will be dominated by special relativistic effects; U can be loosely interpreted 
as the peculiar momentum per particle per mass. If F ^ 1 and f/ > in the wall of 
a void, the wall moves outward with momentum much greater than the gravitational 
attraction inward. We now investigate what happens to this wall. From Eqn ( 2.22| ), the 
"conservation of momentum" equation is 

nw<^ U = -nw (^-^ + j - F ■ (2.29) 

The functions F, n, p, Q, M, w and R' are always positive. If f/ > 0, then at the 
inner edge of the void wall where (p + Q)' > 0, U will be negative — the fluid there is 
decelerated. This deceleration is due to both gravitational and fluid forces. We examine 
the fluid force contribution only. Using Eqns ( p.l2[ ), ( p.6| ) and ( |2.22| ) we find that 



r = -TU^iP±9y, (2.30) 
nwc^R' 

Again, if (p + Q)' > and U > 0, T will be negative and the fluid particles there lose 
their energy per particle mass. 



Suppose the wall has a very large outward momentum so that U/c^ y 2GnM/ (Rc^) > 
1. Then F ~ t//c. Setting p = p/3 and Q = as for a relativistic, non-viscous fluid. 



where we have used the solution for $ from Eqn ( p.l4| ). We consider the deceleration of 



fluid shells at the inner edge of a steep void wall. We can write p'/R' — Apwaii/A-Rwaii 
where Ap^aii is the difference in the energy density over the width of the wall, and A_Rwaii 
is the thickness of the wall. Then, since A/9„aii — Pmax, where pmax is the maximum wall 
energy density, we can approximate F by 



/ PontP 



max 



Ai?wair'. (2.32) 



4 \ / 

Initially, except for the first factor of F^, all factors on the right hand side of the previous 
equation are independent of F. Thus, F oc — F^, which can be a very large damping 
factor! The second factor is proportional to PoutPmax/P^ — (Pout/p)^, so for the mass 
shells on the innermost part of the wall (i.e. those shells with the smallest values of 

p/Pout) this factor can be enormous. The third factor is Ai?„aii~^, so the thinner the 
wall, the faster it will slow down. If the second and third factors change slowly enough 
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with time, then the slow-down time for the wall is roughly independent of its initial value 
of To, since -J-^^dT/T'^ ^ 1 > {p^^t/pfAR^^n'^At/A for Tq > 1. This could be an 
extremely important result, and would imply that the initial peculiar wall velocity could 
never be large enough to cause a void to expand for an arbitrarily long time. However, 
because the second and third factors in Eqn ( |2.32| ) will change with time and depend on 
r, this is only a crude guess. 

As a concluding remark, we note that the wall of a void formed during first-order 
inflation has an enormous outward peculiar velocity. This enormous velocity has been 
thought to cause a void to expand "indefinitely" . However, with such a large deceleration 
of the void wall, it will slow down in a fimte (and possibly small) amount of time. A 
future paper will explore this numerically.^ 

III. Initial Conditions, Boundary Conditions, and Numerical Techniques 

A. Initial Conditions 

The grid used in this code is initially equally spaced: AR{ti) = R{ti)j+i — R{t{)j = 
constant, where the subscript j denotes the spatial grid point number and ra nges from 
j G [0,jB], where Jb is its value at the outer boundary. (Note from Eqns ( p.l9|) and 
( 2.20| ) that in general, Ar ^ constant then). The O*'^ and I'^^-grid points are located at 



R^ -AR{ti)/2 and Ri = AR{t{)/2 respectively. Thus, Rj{ti) = Ri{ti) + {j - l)AR{ti) 
The initial conditions for the functions at t^ are determined as follows. The viscosity 
is set to zero: Q{ti, R) = 0. The energy density p(ti, R) (or the "mass-energy" M(ti, R)), 
the "velocity" U{ti,R) (or V{ti,R)) and the specific internal energy, e(ti,R), are chosen 
as functions of the radius R[ti). (If M is specified initially instead of p, we determine 
p via Eqn. ( |2.6| )). For the cases run in this paper however, we will initially specify the 
fluid to be a polytrope (constant entropy on the initial time slice (II-C)) and e/c^ ^ 1 
or e/c^ -C 1. Using the relations found at the end of Section II-C, e{ti,R) is determined 
once the internal energy e(ti,RB) is specified at the outer boundary {Rb = R{ti)jg): 

e(ti, R) = e(ti, Rb) [ p(ti, R)/p{tu Rb) ] for e/c^ » 1 

e(ti, R) = e(ti, Rb) [ pit,, R)/p{t,, Rb) ] ^^''^ for e/c^ « 1 (3.1) 

We then determine n and p by n(ti,R) = p(ti,R)/{l + e(ti,R)) and p(ti,R) = ( 7 — 
l)n(ti, R)e{ti, R) respectively. Next, we find r (and M if it is not initially specified) by 
integrating outward from r = using the 4*'*-order Runge-Kutta method: 

\ 1/3 

ATmR^dR/T 

l^and M{ti,R) = j\{l + e/c^)R^dR^ , (3.2) 

where we have omitted the (ti,i?)'s for clarity. Finally, $ is found by integrating 
Eqn (|2.26| ) inwards from the outer boundary (again using the 4*'^-order Runge-Kutta 
method) once $(ti, Rb) is specified. 

We are interested in evolving a non-expanding, empty void embedded in a FRW 
universe. We will therefore initially set the energy density outside the void (pout = 
piti^Rs)) to be constant and equal to the spatially fiat FRW value. Then, Hl^^^iti) = 
H'^iti^Rs) = {i/ti)"^ = 8TTGiyc~'^poutiti)/3, where aout(^) = ci{ti){t/t{)^ is the cosmic scale 
factor and cH~^^ is the Hubble radius outside the void. (Note that for the k = FRW 
universe, r = (A7rn(ti))^^^R{ti,r), or a(ti) = {A7rn{ti))~^^^). For p = p/3 and p = 0, 
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^ = 1/2 and ^ = 2/3 respectively. We can quantify the initial size of a void by measuring 
its radius relative to the horizon size outside the void initially: c~^i?waii(^i)/-f^(^tUi) = 
c~^^i?waii(ii)/^i, where -Rwaii(^i) is the void "radius". Following past conventional, we 
loosely equate the Hubble radius with the horizon in the phrase "superhorizon-sized" . 
(Horizon in this context is not to be confused with the particle horizon.) A void is defined 

to be superhorizon-sized if c~^i?waii(^i)/-f^i^t(^i) > 1 ^^(^ subhorizon-sized if c~^i?„aii(^i)/ 



^o"ut(ii) < 1- In addition, Tont{ti) = 1 (see II-B) (or U = ^2GnM/R ) and $out(ti) = 1 
(see II-A). 

In this paper, we consider voids which are initially either compensated or uncompen- 
sated in energy density and which have the following distributions. The "mass-energy" 
function for compensated voids is defined to be 

M{ti,R) = .5pout(ii) [(l + tanhx) + a(l - tanha;)] /2^(ti)/3, (3.3) 

where x = (-R(ti) — -Rwaii('^i))/A-Rwaii(^i), A-Rwaii(^i) is the wall thickness and a is a spec- 
ified constant less than or equal to 1. Because M reaches its spatially flat FRW value 
outside the void, the energy density missing from the void has been put in the wall. For 
uncompensated voids, the energy density is instead initially specified. It is 

p(ti, R) = .5pout(^i)[(l + tanhx) + a(l — tanhx)]. (3.4) 

Here, the energy density missing from the void has not been put into the wall. Therefore, 
the region outside a compensated void will always be a spatially fiat {k = 0) FRW 
universe, whereas the region outside an uncompensated void is a negative spatially curved 
{k < 0) FRW universe. 

The inside of the void is chosen to be homogeneous initially. Since we want it to be 
nonexpanding in the limit that it is empty (p — > 0), we choose the inside of the void to be a 
spatially fiat {k = 0) FRW "mini" universe. We reason as follows. Friedmann's equation 
inside the void is if^ ~ 87rGArc~^pin/3 — kc~'^r'^ / R^ (see II-B), where the subscript 'in' 
denotes quantities inside the void, R = rain and ifm = Um/R = ^^din/ain- In the limit 
that p ^ 0, the inside of the void is non-expanding (din = 0) only if = 0. Since we 
cannot numerically choose pin = (because $in = oo from Eqn( |2.14] )), we would like the 
inside of the void to not expand on time scales that the outside region expands in. This 
is satisfied if pin/Pout ^ 1 because the Hubble time inside the void is much larger than 



that outside the void — H-^^/H^^^ = Jpout/Pm- (As a check, Eqns. ( |2.22| ) and ( p.25| ) 



show that as ?7 — 0, e ^ and p — 0, then U ^ and e — >• 0). Because the void is 
approximately homogeneous, the "mass-energy" inside the void is M{ti, R) ~ c~^pin-R^/3. 



Finally, we set rin(ti) = 1 inside the void since F = — kr'^/c^. The "velocity is then 

U = J2GnM/R ~ c-^rJStvGnp/S. 



We will only consider two types of initial velocity profiles in this paper. The first is 



U = Uq^av = / R (or r(ti,/2) = 1). This specifies that the outward "velocity" 

of each particle is just large enough to compensate for the inward gravitational attraction 
(i.e. the peculiar velocity, t/pEC; is zero). The second is 



U{U, R) = c-^rJs-kGnp/S. (3.5) 



For R < -R„aii(ti) - AR„sii{ti) and R > -Rwaii(ti) + Ai?„aii(ii), F ~ 1. In the wall region, 
however, F > 1. This corresponds to an initial net outward "peculiar momentum" of the 
particles in the wall. 



9 



B. Boundary Conditions 

We set the outer boundary conditions for all times to be those for a homogeneous 
fluid. This specification works well in practice as long as the action is taking place 
away from this boundary. Thus we set ra' = e' = p' = and Q = at j = Jb, where 
Jb is the grid point number for the outermost comoving coordinate. In addition, we 
set = 1, its FRW value (see Eqn (f^.3D. Note that specifying ^jg(t) completely 

eliminates the arbitrariness of the time coordinate, as discussed after Eqn ( p.2|) . The 
present definition sets t to be the FRW cosmic time outside the void. Thus if an initially 
inhomogeneous fluid becomes homogeneous, then from Eqn( |2.26| ), $(t, r) = 1 everywhere, 
and t = constant hypersurfaces correspond to t = constant FRW homogeneous and 
isotropic hypersurfaces. 

It is possible to determine the outer boundary conditions for all t by solving the 
equations with Q = p' = 0. One is then left with two 1**— order ordinary differential 
equations to solve. The boundary conditions determined this way, however, give larger 
errors than the ones shown below, and therefore were not used. We instead integrate 

R = U 

U = - (GnM/R^ + AnGNpR/c^) 

M = -AnpR^U/c^ (3.6) 

using the MacCormack method. We then determine n, e and p from Eqns ( |2.28| ) and 

Eqn (12.181) by setting n^-^ = n^s^i, ej^ = ej^iti) [ nj^^it) /uj^it-,) and Pj^ = {l - 

Finally, reflecting boundary conditions are used at the inner boundary: R^ = —Ri, 
Uq = —Ui, Pq = Pi, ?7.Q = n", eg = e" and Qq = Q", where the superscript n refers to 
values on the n*^ time slice. 

C. Numerical Integration 

The U, e, R, n and M equations are integrated using the 2-step MacCormack 
predictor-corrector method.^a In Appendix B, we give the exact form for the difference 
equations which allows inbound shocks to rebound off the origin. To illustrate the Mac- 
Cormack method, we show the predictor and corrector steps for h as an example. We 
continue to use the convention that n*- is the value of n on the j^^ spatial grid point and 
at the i^^ time step. Suppose we know all quantities on the i^'^ time slice. We would like 
to determine them on the {i + l)*'^ time slice. First we predict the new quantities (with 
forward differencing) using the functional values on the i*^ time slice: 

= n] - At n^m .^^ . (3.7) 

' ' ' ' RfiRj\i-R/) 

After using similar methods to obtain Up^^^, Rp^^-, Mp^^ and Cp^j"^ (and setting Pp"^^ = 

(7 — l)?7,p*."'"^ep*."'"^) for all j, we integrate the $' equation inwards from j = Jb using the 

4*''-order Runge-Kutta method with linear interpolations to determine for all j. We 

integrate again, using the predicted values obtained above (with backward differencing), 
and then average these with the previously predicted values: 
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We obtain the other values in a similar way, and then integrate again to find 

As is well known,lii it is important to choose small enough time steps At to satisfy 
the Courant condition. This condition requires At to be smaller than the time taken for 
sound to cross from any one grid point to the next. The speed of sound for relativistic 

fluids is cs = J {dp/dp)s^^ Using the fact that TdS = d{p/n) +pd{l/n), we find 



OS = J—- (3.9) 

V nw 

As — > oo, w = 1 and we obtain the usual expression for the speed of sound in 
nonrelativistic fluids. Note that for e/c^ ':$> 1, Cg = Ca/7 — 1 = .57c for perfect fluids with 
7 = 4/3. The Courant condition requires the proper time per proper distance to be equal 
to C/cs, where C is the Courant number and is approximately 1/3 to 1/2 for strong 
shocks. Since the inflnitesimal proper time and proper distance are $At and AAr = 
AR/T respectively, the Courant condition becomes At^ = C - R?) / {^i^iics)i )■ 

^ In addition, when Gn ^ 0, regions of the universe expand and contract. We require 
the time steps to be small enough for the functions n, e, and M to change sufficiently 
slowly. We therefore set At" = / {I'i/l^) to be the maximum time step allowed for 
I = n, e and M, where / is a constant less than one. 

After the n*'* corrector step, At^ is calculated for all i, and At""*"^ is set to the smallest 
value obtained: 



(Atr+^ = min( At^ax, C ^'+^ , 7 

V J I max, rf$f(c5)r 



n "'t f * f 



,(3.10) 



J if Giv^Oy 



where Atmax is a specifled upper bound, if desired.E3 

We apply a convergence test to the code for test problems where analytic solutions 
are available. The relative error in q, where q denotes any quantity, is deflned to be 

ei = \qi-q{ri)\/q{ri), (3.11) 

where q{ri) is the exact solution and is the numerical solution. We obtain a global 
measure of the error by deflning 

1 ^ 

^i = T^Eei, (3.12) 

i=l 

where N is the total number of grid points. This error is proportional to the grid spacing 
to some power: Li oc AR^, where s is the convergence rate. If s ~ 2, the code is 
second-order|y-|as desired. These tools have been used previously to test codes in other 
applications .E3 

IV. l^*-Crossing Time for Relativistic Voids 

In this section, we flrst review the standard lore for the evolution of superhorizon- 
sized voids, and then calculate the I'^^-crossing time (the cosmic time taken for a photon 
initially at the inner edge of the void wall to reach the origin) in this picture. We 
then calculate the actual I'^^-crossing time. If at the l***-crossing time the fluid were 
approximately homogeneous and isotropic, distortions from the original void would be 
negligible and a (nearly) FRW universe would result everywhere. 
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It has been suggested that superhorizon-sized voids formed from first-order infla-, 
tion would conformally expand with spacetime during the radiation-dominated period.™ 
Thus, the size of a void at time t would he R = roa(t), where tq is the comoving coor- 
dinate of the void and a{t) is the cosmic scale factor. There are several justifications to 
support this belief. First, small density perturbations conformally expand with space- 
time. Second, vacuum bubbles conformally expand during inflation after an initial short 
growing period. Third, the time taken for the void wall to slow down is expected to be 
enormous because the initial outward momentum of the void wall is enormous. However, 
this reasoning is not enough to conclude that a superhorizon-sized void conformally ex- 
pands with spacetime. First, a void is not a small perturbation in spacetime. Thus linear 
results can not be applied to the description of a void. Second, it is the negative pressure 
that causes vacuum bubbles to expand; a similar configuration having positive pressure 
would instead acclerate inward.^ Third, although part of the wall may still move out, 
the deceleration of the inner void wall is extremely large (see II-D), so that the void can 
still collapse (see VII). 

If a superhorizon-sized void were to conformally expand in spacetime, then the earliest 
time at which thermalization and homogenization can occur is when the horizon is of order 
the size of the void, since this is the expected f^^-crossing time. If the void has comoving 
size ro and the outside Hubble radius (that outside the void) is cifJut(^i)' ^^i^ occurs when 
roa(t) = H~^^{t). For evolution during the radiation-dominated epoch, H~^^{t) = 2t and 
p = p/3 so that the time is of order Ate = t — ti — -f^,^t(^i)(c~^-Rwaii(^i)/-f^i^t(^i))^- If the 
void is much larger than the Hubble radius outside the void (c"^-Rwaii(ii)/ -f^cmt(^i) ^ 1)' 
then the l***-crossing time is very large: At^/ H^^^{ti) ^ 1, independent of the "emptiness" 
of the void. Other authors suggest that the void would fill in with radiation.0 If spacetime 
at the void wall continues to expand as a{t) when the fluid diffuses into the void, the 
comoving radius of the wall is roughly r ^ ro — c dt/a{t) or 

ra{U) ^ i?waii(ti) - c j'^dt'^tjt' (4.1) 

in a radiation-dominated universe. The l*^*-crossing time (i.e. the earliest thermal- 
ization time) then, is when r ~ 0, or when Ate — -5-^,^1 (ti)(c~^i?„aii(^i)/-f^(^t(^i))^) 
which is roughly the same as the time for the horizon to "engulf" a conformally ex- 
panding void.|^ Because R = ra{t), the radius of the void would be given by i? ~ 

{Rwaiiiti) — cH~^^(t{) ^ t / ti )^t/ti, which for nearly all of the time is i? ~ RwaiiiU) the 

void conformally stretches with spacetime. At time t ~ -f^out(^i)(c~^-Rwaii(^i)/-f^oiit(^i))^5 
the void radius decreases quite rapidly to zero. Thus, although the qualitative void evo- 
lution is quite different in these two pictures, the quantitative I'^^-crossing times are not 
because in both the void comoves with spacetime. The important point is that the ear- 
liest possible thermalization time in both pictures (i.e. the l^*-crossing time) is thought 
to be 

Ate ^ H-;,itOic-'R„,nitO/H-;,itOr. (4.2) 



''In the bubble wall, Gjv = 0, \p\' > and p < 0. S Using Eqn(2.22), we see that the acceleration is 
positive, so that the bubble wall moves outward during inflation. But for normal positive pressure with 
p' > 0, the acceleration is negative, so that the wall must accelerate inward. 

•^We thank Michael Turner for this explanation. 
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We will show in this section that the actual I '^^-cr ossing time is remarkably shorter than 
Eqn([4.2|). In doing so we will show that Eqn( [4.1| ) is fundamentally flawed. 

Consider two radially propagating photons A and B. Photon A starts at the inner 
edge of the void wall and propagates inward, and photon B begins at the outer edge of 
the void wall and moves outward. We would like to calculate the distance each travels in 
time At = t — ti, where is the initial time. Using Eqns ( ^.21 ) and the infinitesimal 

coordinate distance traveled by a photon in time dt is dr = cdt^F/R'. Define Pin(^) 
and Pont{t) to be the energy densities inside and outside the evolving wall region of the 
void. (Thus the subscripts "in" and "out" refer only to the undisturbed fluid). Initially, 
Pin(^i) — Pout(^i) = 0' T(ti,R) = 1. Wc also cousldcr non-viscous, relativistic fluids, 
so that p = p/3 and Q{t, r) = 0. 

We will first calculate th e dist ances photons A and B travel in special relativistic 
voids {Gn = 0). From Eqn. (|2.22| ), the fluid acceleration is zero inside and outside the 

void: f/in = f/out = 0. Therefore, R{t,r) = R{ti,r) inside and outside the void, so that 
the infinitesimal distance traveled by photons A and B in time dt is dR = c^Tdt. From 
Eqn. ( p.lO|) , we see that pout and pin are constant in time. Since $out = 1, we find that 



$in(t,r) = $in(ti,r) = {pont{ti)/pin{ti)y^^ from Equ. i^J^. In addition, since T (x p' 
(Eqn. (|2.30|) ), rin(t) = rout(^) = 1- Therefore in time At = t — ti, photon B travels 
outward the distance ARB(t) = Rsit) — -R„aii(ti) given by 

Ai?B = c'I'outroutAt = cAt, (4.3) 

while photon A travels inward the distance ARA(t) = -Rwaii(ii) — RA(t) given by 

Ai?A(t) = C^inFinAt = C [pout (ti) /Pin(ti )] At = C Tout (ti)/rin(ti) At > cAt. (4.4) 

Thus the emptier the void, the farther photon A moves relative to photon B! (It is 
important to note that this is a strictly relativistic effect; for nonrelativistic voids (II-C), 
$(t, r) ~ 1 inside and outside this void so that the distance traveled by photons A and 
B are approximately the same: ARa — cAt = AR^)- Define Ate to be the f ^^-cr ossing 
time (the time taken for photon A to reach the origin). Then from Eqn. (|4.4| ) with 
ARa = -Rwaii(ti), the I'^^-crossing time is 

At, = c-'R^,n{ti) [ Pin(ti)/Pout(ti) . (4.5) 

Since Pout(i^i) > Pm{ti), Ate ranges from Rwaii(ti)/c to zero. Therefore in the limit that 
Piniti) / poutiti) — i> 0, a photon will reach the origin in zero cosmic time. However, calling a 
region a "fiuid" if it is empty (pin = 0) is incorrect. Suffice it to say that the I'^^-crossing 
time can be arbitrarily small. 

We note that if Eulerian, synchronous coordinates were used with the metric ds'^ = 
—c^dT"^ + dR^ + R^dQ^, the distance traveled by photon A or B in time AT would be 
the same: AR = cAT. Thus the time AT, for photon A to reach the origin is AT, = 

R„sii{ti) . (We can also see this using Lagrangian coordinates, since the infinitesimal 
proper time measured by a comoving observer inside the void is dr = ^i^dt = cdT, 
and therefore the time as measured by this observer for photon A to reach the origin 
is Atc = c$inAtc = -Rwaii(ti))- Thus, the quick l^*-crossing time is due to the choice 
of comoving, synchronous coordinates, a choice we do not have for superhorizon-sized 
general relativistic voids {Gn ^ 0) embedded in a FRW universe.^ 

^ If the spatially flat FRW metric is transformed to the Eulerian gauge with synchronous coordinates, 
a coordinate singularity at the Hubble radius (see Appendix A). Thus this gauge and coordinate choice 
cannot be used to describe the evolution of superhorizon-sized general relativistic voids. 
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We now calculate the distance traveled by photons A and B in the same cosmic time 
for a general relativistic void. Again, the infinitesimal coordinate distance a photon 
travels in time dt is dr = c^Tdt/R'. Because the pressure outside the void redshifts due 
to Hubble expansion, $in(t) decreases in time from Eqn( |2.14D (since $out(t) = 1), and 
Ate consequently increases. We first calculate the distance photon B travels. Outside 

the void, rout(t) = 1, and R = ra{t) so that R' = a{t\)\Jt/ti. ii The comoving distance 
photon B has traveled at time t is Ar = ctij a{t\)[^t/t\ — 1), so that the distance photon 

B travels is ARb = Rb - R^A^i) = - l)[i?waii(ti) + 1/2 cH-^t{ts)^-^- 

We now calculate the location of photon A. Assume that the energy density inside 
the void changes negligibly: Pin{t) = constant. (We will address this approximation in a 
moment). Then $in(t) = (poutit) / Pm{U)Y^'^- Because the energy density outside the void 
is redshifted as p oc l/t^(see VI), @ ^m(t) — ^hi{U)\ftiit- In addition, since F oc p' (from 
Eqn. ( p.30| )), rin(t) = rout(^) = 1- Integrating, we find that the cosmic time taken for 
photon A to travel the distance Ai?^ = -RwaiK^^i) 



At = t-ti = c'^AR 



' Pin(^i) 
.Pout(ti) 



1/4 



1 + 



Ra is 



^AR, 



Pm(^i) 

2H-^,{t,) [pontiti] 



1/4 



(4.6) 



and the location of photon A as a function of time is 

Ra = R^At,)-c<!>Uti) J^'^dt'^/tjf 

= /?wall(ti) - C$in(ti)i/ouU^i) 



(4.7) 
(4.8) 



As in the special relativistic case, Ai?^(t) > ARB{t), so that photon A travels farther 
than photon B in the same amount of cosmic time. The f^^-crossing time relative to the 
initial outside Hubble time is then {ARa = -Rwaii(^i)) 



Atr 



C ^i?wall(ti) / Pin(i^i) 



1/4 



Hontiti) \Pont{ti) 



1 



^-Rwall(ii) / Pin('^i) 



1/4- 



2-ffout(^i) 



,Pout(ti) 



(4.9) 



In addition, if the more stringent condition c ^-Rwaii(ii)/-f^out(^i) < '^'in(^i) holds, then 

Atr , C"^i?wall(ti) ^ /Pout(ii) V''^ 



< 1 



when 



< 



Pin(ti) 



(4.10) 



— the minimum thermalization time (i.e. the P*-crossing time) is less than the initial 
Hubbl e ti me ou tsid e the void! (Note that the general and special relativistic results 
(Eqns (^^ ) and (|4.5| )) are equal in this case, since the energy density outside the void is 
constant during this time). 

Before discussing further implications, we find the condition for which Eq n(|4.6| ) is sat- 
isfied; it was derived under the assumption that p-mit) — Pm(ti). Using Eqn (|2.10|) , inside 

the void \p/p\ = A{R^R)' /{3R^R') = AR/R (since R = ra{t)) so that if pi^ ^ constant, 

then R ~ constant (in time) inside the void. Because U"^ = (i?/$)^ = 2GnM/R ~ 
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c ^(STrGArpin/S)/?^ inside the void, the fractional change in the radius R over time scale 
At is roughly 



^ Ai? R^t 



R R \ Pout(t) 




^^^^-i/out(t)$in(t)At ^ K'{t-Ml"i 77^' (4-11) 



where we have used the fact that Hj^{t{) = 2ti, y^Pin(t)/Pout(i)-f^out(i) - '^'in^(ii)-^^out(ii) 
and = (Pout(ti)/Pin(ti))'/'. We require / < 1. W riting t = At + Eqn(pn|) 

becomes At = f^l{U)H^,^{k)[l + ^1 + l/(/$in(ti))2]. Since <l>in(ti) > 1, we have 
/<I>in(ti) > 1. The condition for which the density inside of the void remains approxi- 
mately constant during time At then, is At/H~^^{ti) < 2p^f^{ti) ~ ^f^{ti). We now com- 
bine this with Eqn( |4.6| ) to find the maximum allowed void size ARa/ H~^^{t\) given $in(ti). 
We find c-^ARA/H^iit,) < .5<l>i„(ti)[-l + y^l + 4(/$i,(tO)2] ~ f^M- Eqn (pj) is 
then valid when /\Ra/ H'^^it^) < f^fjt^), and Eqn. (|0| ) is valid when 

C-'i?wall(ti) / H~^,{t-;) < \/pout(ti) / Pm(ti). (4.12) 



Note that Eqn( [4.10D is automatically satisfied. 



The quick l^*-crossing time might seem completely counterintuitive. How can a pho- 
ton travel a distance much larger than the Hubble radius in less than a Hubble time? The 
answer lies in describing how one measures the size of an object which is not a small per- 
turbation in spacetime. If size is measured circumferentially, then the void is enormous 
because its circumferential size is 27ri?„aii(ii) ^ -f^out(^i)- (If spacetime were static, then 
it would take a photon time At ~ 27ri?„aii(^i)/c to encircle the void). However, if size 
is measured radially (the time taken for a photon to cross the object if spacetime were 
static), then the void is measured to be very small. If fact, an interesting comparison 
can be made to measuring the size of a black hole, an overdense region. If one measures 
its circumferential size, it is small (or at least finite), but its radial size is infinite. 

The l***-crossing time for a superhorizon-sized void was previous calculated incorrectly 
because t was assumed to be the proper time outside and ins ide the vo id; i n our notation, 
it was implicitly assumed that $(t, r) = 1. Comparing Eqns( ^nD and (|4.7|) , we indeed see 
that the factor $in(ti) > 1 is missing from Eqn (]4.1| ). Note in addition that the factor of 

^Jti/t in Eqn( [4.7|) does not come from spacetime expanding at the wall, but rather from 

the density outside the void redshifting, causing $in(t) to decrease. In fact, spacetime 

is roughly non-expanding at the inner wall edge. If c~^-Rwaii(^i)/-f^<^t(^i) > "^111(^1) ^ 1, 
the actual position of photon A is approximately constant in time until the last moment: 
Ra — -Rwaii(^i) until t ~ ti + Ate (see Eqn( [4.8D ). Thus, spacetime at the inner edge of the 
void wall is neither expanding nor contracting. 

We note the interesting fact that if a photon inside the void is in thermal equilibrium 
with average frequency z/, its frequency is Tin. If it moves outside the void, its frequency is 
blue-shifted to u' = z/$in = uiTont/Tin) = Tout, the average frequency of thermal photons 
outside the void. Thus, this photon is automatically in thermal equilibrium outside the 
void, so that an outside observer could not detect the void's initial presence unless non- 
thermal photons came out.^ This is essentially because the fluid is relativistic, or that p, 

^Because entropy is created at the shock (see VII), the photon with j/Tin would actuall have a shghtly 
higher frequency than thermal outside the void. 
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p, n and e depend only on the temperature (II-C). We can understand why the photon 
is blue-shifted by way of comparison to an (overdense) black hole. Suppose an observer 
falls into a black hole ticking off photons at a fixed frequency. As this observer crosses 
the event horizon, the frequency of the photon emitted lasi gets redshifted to infinity as 
observed by a stationary observer outside the black hole.t3 The opposite effect happens 
for a photon emitted from an underdense region. Because a photon leaving a void enters a 
region with a much larger gravitational potential, the frequency instead gets blue-shifted. 

In conclusion, the f^^-crossing time for a super horizon- sized relatiyistic void embed- 
ded in a FRW expanding universe is given by Eqn (if Eqn('gl2|) is satisfied), and 



depends sensitively on the quantity -Rwaii(^i)/-f^out('^i) {Pm{ti) I PovLt{ti)Y^^ ■ We emphasize 
the important point that if c"^i?waii(ii)/-f^out(^i) < (Pout(ii)/Pm(ii))^/^ = T^ntiU) / TiniU) , 
then the I'^^-crossing time is less than the outside Hubble time: /S.tc/ H^^^{ti) < 1. 

V. Pressureless Non- Viscous Voids 
For the special case when the pressure and viscosity are zero, the equatiqns of motion 
can be solved analytically. These give the Tolman-Bondi dust solutions,u'Ej which we 
will review here briefly. It is important to study the pressureless case not only as a test 
problem, but also to see how removing fluid forces affects the evolution of a void. (Because 

p = Q = , the particles move only under gravitational forces: $^^[7 = —GnM/B? (see 
Eqn iV^m )). 



Since $ = (Eqn ( p.26| )), we set ^(t, r) = 1. The mass contained within r re- 
mains constant: M(r) = A-rmR^dR/T , since M = (see Eqn ( |2.9|)). And from 
Eqn(p30|), r(t,r) = r(ti,r). For R' ^ 0, p = n can be found from Eqn(|2.19|): p(t, r) = 
piti,r) Rj tj, r)' R{ti, r)"^ /\R{t, r)'R{t, r)^]. The generalized pressureless Friedmann equa- 
tion, Eq. (|2.12|) , now becomes 



R' = c 



Fir 



i2 



+ 2GNM{r)/R. (5.1) 



The quantity — GnM/R = c^(r^ — l)/2 is conserved during evolution and can be 
interpreted as the generalized total energy. Eqn ( p.l[ ) can be integrated for a shell of 
radius r: 



R = Y2_i (coshr^-l), t = To{r) + j^^^—j^ [smhr] - r]) , for r(r) >1 



R = (9GjvM/2)^/^ (t - To{r) f^ , for r(r) = 1. (5.2) 

These are the Tolman-Bondi solutions. 

In these models, shell-crossing can occur. This happens when two adjacent shells 
(labeled by r and r + dr) occupy the same position so that R' = and p — >■ oo. This may 
lead to a non-unique continuation of the solution, BZl and thus computations have to be 
stopped. This problem is beheved to occur because the pressure has been artificially set 
to zero. It is generally thought that adding pressure would prevent this situation from 
occurring .0 As will be seen, when r(r) > 1 in the void wall region, shell-crossing does 
occur. The addition of enough artificial viscosity can prevent this from happening, how- 



ever. Even though the viscosity given by Eqn ( 2.27 ) was designed to stabilize numerical 



shocks, it has been found to prevent shell-crossing. 
''The author thanks T. Piran for this information. 
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Table 1: Convergence test for Tolman-Bondi model 



AR{ti) 


Li 


s 


Lla 


s 


Lib 


s 


8 


1.2 X 10-3 




8.33 X 10-^ 




8.55 X 10"^ 




4 


4.13 X 10-^ 


1.48 


2.26 X 10-^ 


1.89 


2.16 X 10-^ 


2.0 


2 


1.55 X 10"^ 


1.32 


6.09 X 10-5 


1.82 


5.35 X 10-5 


1.97 


1 


6.56 X 10-5 


1.24 


1.81 X 10-5 


1.75 


1.41 X 10-5 


1.92 



For the numerical simulations in this section, we set Gat = 1, c = 1, ti = 1, p{ti, R) = 
0, eiti^Rs) = 0, C = .3, / = .005 and 7 = 5/3. Thus, the initial energy density and 
Hubble radius outside the void are 47rpout(^i) = 2/3 and H~^^{t-^ = 3/2, respectively (see 
III-A). 

We first examine the situation in which each mass shell's velocity initially compensates 



for the gravitational attraction inwards: T{t[,R) = 1. Then U{ti,R) = \J2GnM / R. 
Figure 1 shows the energy density versus R R{ti, Rb)/ R{t, Rb) (= R Rjg{t{)/ Rj^it)) for 
a superhorizon-sized compensated void with c"^i?waii(^i)/-f^(^t(^i) = 333, Ai?waii(^i) = 15, 
a = .001, AR{ti) = 2.5 and fc^ = 0. We show the analytic and numerical results at 
times t = 1, 10, 100, and 300 where Rjg{t) = 1002, 4651, 2.16 x 10^, and 4.49 x 10^ 
respectively. The triangles and squares are the numerical and Tolman-Bondi solutions, 
respectively, although they are difficult to distinguish because the numerical results agree 
so well with the analytic results. By t ~ 300, the density everywhere is approximately 
constant; there is hardly a trace of the void's initial presence. Identical results are 
obtained for subhorizon- sized voids. In addition, an initially uncompensated void evolves 
similarly. 

In Table 1, we show the results of a convergence test for g = p applied to the same 
initial conditions as in Figure 1, but for variable AR{ti) (see III-C). (We only do this 
test for p, because the accumulated error in M and R are much smaller). We set up 
analytic conditions initially and integrate until t = 1.5. Because the inner grid point 
or two ends up being the numerical culprit for non-second order convergence, we also 
calculate Lia = Hf=2, ^-nd Ln, = j^^jq Y^iLio ^i, where Cj is the relative error. (This 
is because the code consistently underestimated p at the innermost few grid points). To 
the right of each global error estimate, the convergence rate is shown (Lj oc AR{tj)^). 
We see that the convergence rate for Li is less than second-order, whereas that for Li is 
nearly second order. 

We have just seen that if f/ = f/cRAv (r{U,R) = 1), the void disappears. What 
happens when U > ?7grav (r(^i5 -R) > 1) in the wall region? Recall that this corresponds 
to a net outward peculiar velocity (II-D). Figure 2 shows the result for a compensated 
void with c-^R^^n{ti)/H-^t{ti) = 333 , AR^^nik) = 15, a = .001, AR{ti) = 2.5, and 

= 0. We choose the initial velocity to be given by Eqn(|3.5|). Therefore, T{ti,R) = 1 
everywhere except in the wall region, where r{ti,R) > 1. ^Again, the triangles and 
squares represent the numerical and analytic solutions, respectively. The initial time 
[ti) and t = 1.02, 1.048 are shown. Again, the difference between the numerical and 
analytical results are small except at t = 1.048, where shell crossing occurs in both 
solutions, a good check on the code. The comoving radius of the shell with the highest 
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density, rsheii(^), remains approximately constant in time. Because each shell in the wall 
has constant total energy c^(r^ — l)/2, the wall expands outward. (Since U > 0, R 
increases. But since the total energy is constant and M(r) is constant, tf must increase). 
Identical results are obtained for subhorizon voids. 

Figures 3a and 3b show the density as a function of R R{ti, Rb) / Rit, Rb) for a 
subhorizon-sized, shallow, uncompensated and compensated void, respectively, at t = 
1, 50, 100, 400 and 1000. Note that shell-crossing would have occurred at t = 116 and 2.1 
for the uncompensated and compensated voids, respectively. Here, c~^Rwa.n(ti)/H~^^(ti) = 
.0067, Ai?„aii(ti) = -001, AR{ti) = .0001, a = .5, and the velocity U{ti,R) is given by 
Eqn (U). In addition, for Figure 3a, fc^ = 4 and Rj^{t) = .035, .49, .79, 2.13 and 4.26, 
while for Figure 3b, k"^ = 8 and Rjg(t) = .035, .48, .76, 1.9 and 3.5. As can be seen, the 
initial perturbation grows with time and eventually forms a thin, dense shell. Again, the 
comoving coordinate for the shell with the highest density is approximately constant in 
time after the shell has formed. As the shell travels outward, it pushes mass in front of 
it, producing a shock. This situation is similar to that of a fast car colliding with slower 
ones; although the faster car is not allowed to move through the slower cars, momentum 
is still transferred to them. Note that the initially compensated perturbation forms a 
thick shell more quickly than the uncompensated perturbation, although it then proceeds 
to grow more slowly. 

As long as the voids formed remain, subhorizon-sized, they will eventually grow ac- 
cording to a known similarity solution.c2 An initially compensated (uncompensated) per- 
turbation in an expanding FRW p = universe will eventually form a dense, thin shell 
that expands outward as -Rsheii = -R(^, ''^sheii) t^^^ (t^^^). In Figure 3c, we show the 
position of the void wall versus time for the results of Figures 3a and 3b. The trian- 
gles and squares connected by lines are the numerical solutions for the compensated and 
uncompensated cases, respectively, and the dashed and dotted lines are the self-similar 
solutions for the compensated and uncompensated cases, respectively. As can be seen, the 
initially compensated perturbation approaches the similarity solution more quickly than 
the initially uncompensated perturbation, but for t > 800, both solutions are self-similar. 

VI. FRW Homogeneous Cosmologies 

We now test our code against the exact FRW homogeneous and isotropic solution 
for relativistic fluids with p = p/3. As discussed in II-A,B and III-A, the FRW solution 
is R{t,r) = rait) = R{t;,r)(t/t;f where ^ = 1/2. In addition, AT:GNp{t) = 3cV(8t^) 
and U = R/{2t). For the numerical results shown in this section, we set = 1, 
c = 1, 7 = 4/3, C = .3, ti = 1, e{tuRB) = 10^ and c-^RBit-,)/ R-^it-,) = 250 so that 
47rp(ti, Rb) = 3/8. 

Figure 4 shows the relative error in p (III-C) for a simulation with AR(ti) = .5 and 
= 0. The analytical solution was set up initially, and the code was run until t = 1.1. 
The solid, dotted and dashed lines are for / = .01, / = .005 and / = .0025. The relative 
error at the outer boundary is seen to be very sensitive to /; it is .05%, .15% and 2% for 
/ = .0025, .005 and .01 respectively. Note also the underprediction of p at the innermost 
few grid points. 

Table 2 shows the accumulated error at time t = 1.1 for simulations with fc^ = 4. [| 
Again, Lib = jv^Yq Si^io where is the relative error for q = p. In the first 7 columns, 
we show the resuhs for = 25, 50, 100, 200, 400, 800 and 1600 (i.e. AR{ti) = N/250) 

*When evolving voids, viscosity is absolutely necessary. It is therefore important to see how it affects 
the solution where the fluid is approximately homogeneous and isotropic. It turns out that it is virtually 
unaffected by Q 7^ 0, as it should be. 
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Table 2: Ultrarelativistic Homogeneous convergence test 





Number of grid points (N) 




Li 






25 


50 


100 


200 


400 


800 


1600 


/ 


AR^O 


s 




.016 


.010 


.0073 


.0058 


.0050 


.0046 




.01 


.004 




Lib 


.0071 


.0056 


.0049 


.0045 


.0044 


.0043 










Li 


.012 


.0067 


.0038 


.0023 


.0016 


.0012 




.005 


.00085 


2.00 


Lib 


.0015 


.0012 


.0010 


.00091 


.00088 


.00087 










Li 


.011 


.0059 


.0031 


.0017 


.00096 


.00060 


.00042 


.00025 


.00025 


1.77 


Lib 


.00057 


.00041 


.00032 


.00027 


.00026 


.00025 


.00025 









and for / = .01, .005, and .0025. For a given value of /, as increases, Li and Lib 
approach the same constant value even though Li starts out much larger. After this 
constant value has been reached, Li and Lib remain unchanged when AR{ti) is further 
decreased, even though the relative error near the origin improves. This is because 
beyond a certain point, all the accumulated error comes from the outer boundary, which 
is immune to changes in Ai?„aii(^i) (see III-B). 

Knowing that a given value of / limits convergence of the code, we can calculate the 
convergence as a function of the asymptotic value of Li. We assume that Li{AR(ti) 

0) oc 7 . Column 9 gives the estimated value for Li{AR{t[) —>■ 0), and column 10 

estimates the value of s. We see that s is nearly 2, which means that convergence in / is 
nearly second order given a small enough value for AR{ti). 

VII. Numerical Evolution of Voids 

A. Nonrelativistic Fluids 

In this subsection we examine the evolution of voids composed of nonrelativistic par- 
ticles in zero gravity {Gn = 0). If T is the fluid temperature and /i is a fluid particle's 
mass, T/n <^ 1. Since oc Gn~'^, these voids are subhorizon-sized. For the simu- 
lations in this subsection, we set = 0, c = 10^°, = 1, C = .3, 7 = 5/3, C = .3, 
Anp{ti,RB) = 2/3 and e{ti, Rb) = 1. 

We start with the shock tube problem, a standard test of 1-D slab codes.ii In addition, 
it provides insight into the dynamics of collapsing voids. In a shock tube, the fluid is 
initially at rest and is separated into two regions with different pressures and densities. 
The pressure discontinuity produces a shock wave which propagates into the low pressure 
region and a rarefaction wave which propagates into the high pressure region. An analytic 
similarity solution exits for a perfect fluid with slab geometry. It does not exist for the 
spherically symmetric geometry however.0 Far from the origin however, the sDherically 
symmetric solution approaches the slab solution for small times and distances.Eil We will 
thus set up a spherically symmetric shock tube by evolving an uncompensated void far 
from the origin, and compare the results to the exact slab similarity solution (briefly 
reviewed in Appendix D). 
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In Figure 5 we show the resuhs for the shock tube problem with U{ti, R) = (or 
T{ti,R) = 1), i?waii(ti) = 20, AR^^niti) = -01, AR{ti) = .01, a = .01 and P = 5. (We 
do not take the pressure gradient to be discontinuous initially, because the solution in 
the original wall area is not as accurate then.0) We plot the pressure, number density, 
velocity and specific energy as a function of the position R at ti and at t = 1.15. The 
triangles connected by lines is the numerical solution, and the dashed lines is the slab 
similarity solution. At t = 1.15, a strong shock wave (located at i? ~ 19.68) moves 
inward and a rarefaction wave (between 19.85 < R < 20.14) moves outward. Note that 
the shock is spread out over /c^ ~ 4— 5 grid points. The numerical and analytical solutions 
are seen to agree well in this limit, because the spherical geometrical effects are small 
((-RshocA:(^) — -Rwaii(^i))/-Rwaii(^i) ~ .4/20 = .02). The distortiou of the vclocity distribution 
is due to the small geometrical effect; the shock gets slightly stronger directly behind the 
shock due to the smaller effective volume AnR'^AR those mass shells occupy relative to 
shells further back. An important point to emphasize is that although initially the fluid 
is everywhere stationary {U{ti,R) = 0), it acquires a net momentum to the left in the 
wall region. 

In Figure 6a we show the long-term results for an initial conflguration with -Rwaii(^i) = 
1 but otherwise identical to Figure 5. The pressure, number density, velocity and speciflc 
energy are shown at the initial time ti = 1 with dashed lines, and long after the collision 
at t = 2.5 with triangles and connecting lines. Although it is not shown for clarity, 
the fluid conflguration before the shock rebounds at the origin is similar to Figure 5: a 
shock heads toward the origin and a rarefaction wave moves away from the origin. When 
the (spherical) shock crashes into the origin, the volume effect proves harsh as the fluid 
collides with itself in a vanishingly small volume. This causes the pressure at the origin 
to become very large in order to repel the fluid (not pictured here). Note that the fluid 
at the shock as well as far behind it must be repelled, since it is all moving toward the 
origin.^ When the dust settles, we flnd that a weak shock has rebounded back. This 
outward-moving shock can be seen at i?(t=2.5) ~ 1.3. Due to the volume effect however, 
this shock will become weaker as it moves outward further. Also seen at t = 2.5 is the 
original (outgoing) rarefaction wave located between 1.3 < i? < 2.4. Note that at t = 2.5, 
the fluid is (and will approximately remain) at rest near the origin because the velocity 
is zero and the pressure is constant. However, a large distortion has been left behind in 
the fluid in the form of low density and high kinetic energy. This consists entirely of the 
fluid originally in the void. The low-n, high-e values for the flrst 5-6 grid points, however, 
is artiflcial. This is a consequence of using VonNeumann-type artiflcial viscosity called 
"wall-heating" , and is caused by the collision between 2 shocks.^ 

In Figure 6b we show the initial and long-term pressure, number density, velocity and 
speciflc energy as a function of the radius for a compensated, nonrelativistic void. The 

■'This is primarily because we calculate n via h oc {R^dR)~^ rather than n oc [r'^dr)^^ 

'^It is this reversal that requires a very robust difference scheme. All of the "ob vious" difference 

schemes failed after the inbound shocks reached the origin. Setting / = in Eqn (|2!20|) and using those 

difference schemes listed in Appendix B are the very necessary requirements. 
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Table 3: Nonrelativistic collapse times 



Void type 


Ai?wall(ti) = .02 


Ai?„all(ti) = .04 


Ai?wall(ti) = .08 


A-Rwall(^i) = .12 


Collapse Time Ate 


Percent Change: (Ate — Ate 


;.02),„)/Ate(.02),n 


uncompensated 


.38 


.41 


II 8% 


.45 1 


1 18% 


.50 1 


1 32% 


compensated 


.15 - 61% 


.20 


-47% 


.27 


-29% 


.32 


-16% 



initial distribution is identical to Figure 6a except for the energy density, and Ai?waii(^i) = 
.02. Again, the fluid configuration at tj = 1 is shown as dashed lines and that at t — 2.5 
is shown as triangles connected by lines. An inbound shock is again formed from the 
initial pressure gradient at the inner edge of the void wall. Unlike the uncompensated 
case however, a weak outgoing shock wave is formed instead of a rarefaction wave. At 
t — 2.5, it is located at i? = 3.2. Like the uncompensated void, the pressure at the origin 
becomes very large after the shock collides with itself there, and a weak shock rebounds. 
At t = 2.5, this rebounded shock is located at i? ~ 2.3, which is farther out than that for 
the uncompensated void {R ~ 1.3). This is because the inbound shock produced from 
the compensated case is much stronger than for the uncompensated case. And again 
there is a large distortion left near the origin containing all the fiuid initially in the void, 
although it is somewhat different spatially. 

We now compare the collapse times for uncompensated and compensated voids of 
varying wall thicknesses. (The collapse time is defined to be the time taken for the shock 
to reach the origin). We set -Rwaii(^i) = 1, a = .01, /c^ = 3 and AR{ti) = .01. Table 3 
shows the results for Ai?waii(^i) — -02, .04, .08 and .12, where we calculate the percent 
change by comparing the collapse time with the uncompensated Ai?waii(^i) = -02 collapse 
time of At = .38. As Ai?„aii(^i) increases, the collapse time increases. In addition, the 
collapse time for compensated voids is substantially smaller than for uncompensated 
voids. 

In conclusion, nonrelativistic voids with zero gravity collapse in the form of a shock, 
the strength of which depends on the details of the void wall. Some time after collapsing, 
the fiuid is virtually at rest everywhere, with n and e inhomogeneous near the origin. 

B. Special Relativistic Fluids 

In this subsection we consider the evolution of special relativistic voids with T//x ^ 1, 
where T is the temperature and n is the mass of a fiuid particle. Thus, e/c^ ^ 1, and we 
set = 0, c = 1, ti = 1, C = .3, 7 = 4/3, e(ti, Rb) = 10^, U{ti, R) ^ (or r{ti, R) = 1) 
and 47rp{ti, Rb) = 3/8 in this subsection. 

In Figure 7a, we show M, F, 47rp and $ as a function of R for a relativistic shock tube 
problem. Here i?waii(ti) = 1, = 3, Ai?„aii(ti) = -02, AR{ti) = .01 and a = 10"^. As in 
the nonrelativistic case, an inbound shock is formed. (This is observed most easily in the 
plot of Attp versus it!). However, an outgoing "rarefaction wave" is not observed, as it is in 



21 



the nonrelativistic case (see Figure 5). If a wave were to propagate outward, it could go 
at most the distance a photon would travel. But from Eqn( [4.3|) , we know that a photon 
starting from the outer wall edge only moves the distance ARb = .04 (.065) in time 
At = .04 (.065). Since this is of order the grid point thickness, if an outbound wave were 
present, it would not be observed at this time anyway. On the other hand, an inbound 
photon starting from the inner wall edge would travel the distance ARa = .4 (.65) (from 
Eqn ([4.4|) ) in time At = .04 (.065), since Pout(^i)/Pm(^i) — 10^. Reexamining Figure 7a, 
we now notice an important result; the inbound shock's position is approximately equal 
to photon A's location — the shock moves inward at roughly the speed of light. This is 
actually not so surprising, because the speed of sound for a perfect fluid with p = p/3 is 
.57c (III-C). 

Consider next voids compensated in energy density. Because of relativistic-particle 
diffusion, we expect the inbound shock to again travel at approximately the speed of light. 
Since the value of $in(ti) does not depend on the functional form of p in the void wall 
(Eqn ( p.l4|) ), the shock should move the same distance per time as for the uncompensated 
void. We ran a compensated void simulation with the same initial conditions as from 
Figure 7a (except for p in the void wall), and found that at t = 1.04 and 1.065, the 
compensated shock is ahead by only AR = .05. However, this can be accounted for by 
the slightly different values of p initially at the inner edge of the wall. Thus, special 
relativistic compensated and uncompensated voids collapse at approximately the same 
speed, unlike nonrelativistic voids (see Table 3). 

In Figure 7b, we show the numerical results for an uncompensated void with -Rwaii(^i) = 
1, p = 3, a = 10"^ and AR{ti) = .01, and for Ai?^aii(ti) = -02, .04, .06 and .1. It is 
clear from this figure that in all four cases the shock reaches the origin at Ate ~ -08; 
Ate is approximately independent of Ai?waii(^i)- (Compare this with Table 3). Using 
Eqn (|4^), we estimate the time for light to reach the origin at Ate — Rwaiiiti)/^ — -09, 
in agreement with the simulations.! (During this time, photon B would only move out- 
ward the distance ARb = Ate — -09). The value of A_Rwaii(^i) however, does influence 
shock formation-time and strength. As Ai?waii(^i) increases, the shock formation time 
increases and the shock strength decreases. 

We now examine the what happens to the void after the shock collides at the origin. 
In Figure 7c, we show the numerical results for an uncompensated void with -Rwaii(^i) = 1, 
Ai?waii(ii) = -02, a = 10-^ fc2 = 8 and AR{ti) = .01. □ We show p, n, U and e as 
a function of R at the initial time t[ and at times 1.04 and 3.0, where the second and 
third times are before and after collision at the origin, respectively. At t = 3, p' ^ and 
f/' ~ near the origin; the fluid there is roughly at rest. Two weak outgoing waves are 

'Because "I>in actually increases slightly due to the dissipation of energy at the shock, the actual travel 
time for the shock to reach the origin, Ate, is decreased slightly. 

'"After colliding, a very weak shock rebounds. The generalized functional form for the artificial 
viscosity, however, was derived in the strong shock limit. Thus, a large value of k"^ was needed to 
maintain numerical stability after the collision. A new-functional form for Q will have to be used in 
future simulations to stabilize the weak outgoing shock.Lj 
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observed: the shock (at i? ~ 1) and a "rarefaction wave" (between 1.1 < < 2.1). In 
addition, in contrast with nonrelativistic voids, n' ~ and e' ~ near the origin.^ This 
result is expected; since n oc p^^^ oc and e oc p^^"^ oc T (from the discussion following 
Eqns ( p.28| )) for a non-viscous fluid, if p' ~ 0, then it follows that n' ~ and e' ~ 0. 
This is a consequence of the fact that p, p, n and e depend only on the temperature for 
e/c^ ^ 1 and Q = 0. We therefore find that after the collapse, there is only a small trace 
of the void's initial presence! 

We conclude that a special relativistic void collapses in the form of a shock which 
travels at approximately the speed of light into the void. Thus, the collapse time is of 
order the l^*-crossing time. Some time after the collapse, the fluid becomes approximately 
homogeneous and isotropic everywhere. 

C. General Relativistic Fluids 

In this subsection we study the evolution of general relativistic voids for T/yU 3> 1, 
where T is the temperature and fi is the mass of a fluid particle. We set = 1, c = 1, 
ti = 1, C = .3 and 7 = 4/3. Therefore, the outside Hubble radius and density (that 
outside the void) are cif^t(ti) = 2 and AtiG ^ p{t\^ Rb) = 3/8. 

In Figure 8, we show the pressure, number density, velocity and specific energy 
for a general relativistic void at the initial time ti and for t = 1.04 and t = 1.065. 
The initial conditions are identical to those of Figure 7a, except that k"^ = 4 and 
U(ti,R) = f/cRAv = ^J^GnM/ R (or r{ti,R) = 1). Because the void is subhorizon- 
sized {c~^Rwiin(ti)/H~^^(ti) = 1/2), its evolution looks virtually the same as that for 
the special relativistic void shown in Figure 7a; at the void wall, the gravitational 
force is M/R^ + AnpR/c^ ~ SnpR/S ~ 10~^, which is much less than the fluid force 
T^p'/{4pR') ~ (4Ai?(ti))-i ~ 12 (see Eqn. (p])). 

In Figure 9a, we show the pressure as a function of R R{ti, Rb) / Rit, Rb) (= R{t,r) 
Rjai^'^ / Rjai^)) ^ general relativistic superhorizon-sized void at ti, t = 2.0 and t = 8.0 
with r(ti,i?) = 1, i?waii(ti) = 50, Ai?^aii(ti) = 1., A;2 = 4, a = 10-^ e(ti,i?B) = 10^ and 
AR(ti) = .5. (This void is 50/2 = 25 times the outside Hubble radius). In addition, 
Rjg{t) = 100.2, 139.1, and 271.8. Even though G^ 7^ here, a strong inward shock still 
forms. This is because particles are diffusing into the void, having been accelerated away 
from the high-pressure wall. (Note that at the wall, the acceleration due to gravity is 
GnM/R^ + A'npR/c^ ~ 10, whereas that due to the fluid force is only T'^p' /{ApR') ~ .25. 
This small relative amount however, is enough to form the shock). Note that because 
of expansion, the pressure outside (and to a lessor extent inside) the void redshifts. By 
t = 8, the pressure outside the void has redshifted from p ~ 10~^ to 2 x 10~^, the expected 
amount since p oc 1/t^ so that Pout(8) — 10~^/64 ~ 2 x 10~^. At the same time, the 
pressure on the inside has only redshifted from p ~ 10^^ to 5 x 10^®, because the inside 
Hubble time is larger than that outside: {ti) / H~^^{ti) = \J pontiii) / Pmiti) = 100. 
From Eqn(^), the P*-crossing time is At^ ~ 50/10(1 + 25/(2 x 10)) = 11.25, in rough 

"For the innermost 7-8 grid points, e and n are too large and too small, respectively. This is again 
due to "wall- heating" . 
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agreement with the numerical collapse time of Ate — 8. Thus, the collapse time is 
found to be approximately equal to the l^*-crossing time — the the shock moves inward 
at roughly the speed of light. This is not surprising however, because the speed of sound 
is .57c (III-C). 

Figures 9b and 9c show the pressure as a function of R R(ti, Rb) / R(t, Rb) for a 
general relativistic superhorizon-sized void with c^^-Rwaii(i^i)/-f^(^t(^i) = 25, r(ti,_R) = 1, 
k"^ = 4, Ai?„aii(ii) = 1, e{ti,RB) = 10^ and AR{ti) = .5. In addition, a = 10~^ 
and Rj^{t) = 100.2, 125.1, and 148.7 for Figure 9b, and a = 10-^° and Rj^it) = 
100.2, 102.6, and 104.4 for Figure 9c. The numerical collapse times in Figures 9b and 9c 
are Ate = 1-3 and .09, respectively, which are both smaller than the outside Hubble time. 
(Note that Ate = -09 is 1/20*^^ the outside Hubble time). Since equals 31.6 and 

316, respectively, the l*^*-crossing times from Eqn( ^4.9| ) are Ate — 50/31.6(1 + 25/31.6) ~ 
2.2 and 50/316(1 + 25/316) = .17, respectively. The I'^^-crossing times are larger than 
the numerical collapse times because $in(t) increases during the collapse due to the 
dissipation of energy at the shock. Thus, the gain in entropy at the shock only makes 
the collapse time shorter. For example, in Figure 9c, $in(t) increases from its original 
value of $in = 316 to $in ~ 450 — 500 during the collapse. Using = 500, we would 
predict Ate ~ -l, which is roughly correct. 

We note however, that in Figure 9c (and to a lesser extent in Figure 9b), only part of 
the void has been filled in at the collapse time te = ti + Ate- Thus, thermalization and 
homogenization has not been achieved by te- We note from Figures 9a,b and c that as 
the initial relative energy density inside the void decreases, the fraction of the void filled 
in at the collapse time decreases. However, since the energy density inside the somewhat 
filled void is still much less than Pout(^c); ^(^c) inside the "void" is still much greater 
than one. And because T{te) is also larger than one inside the "void", the distance light 
can travel is still greater inside than outside the void. Thus, although the void has not 
yet homogenized, this may only take an additional small amount of time. 

Figures 10a and 10b show the pressure versus R R(ti, Rb) / R(t, Rb) for a. superhorizon- 
sized uncompensated and compensated void, respectively. Here, T(ti,R) = 1, -Rwaii(^i) = 
500, Ai?„aii(ti) = 10., A;2 = 4, e{ti, Rb) = 10^ AR{ti) = 5., and a = 10"^°. (These voids 
are 250 times the outside Hubble radius). The pressure is shown at ti = 1 and t = 1.7 
for each void, and at t = 2.1 and t = 2.0 for the uncompensated and compensated voids, 
respectively. In addition, in Figure 10a, Rjg{t) = 1002, 1288, and 1424, while in Figure 
10b, Rjg{t) = 1002, 1307, and 1418. The important point to note is that the shock 
reaches the origin at t ~ 2.1 for both voids. This is due to the diffusion of particles 
into the void, and does not depend on the compensatedness of the void because $in(ti) 
depends only on Pin(ti)/Pout(^i)- (This is similar to that for special relativistic relativistic 
voids (VII-B)). This roughly agrees with the predicted collapse time Ate = 500/316(1 + 
250/(2 X 316)) ~ 2.2 from Eqn(|]9]). Note also that at t ~ 2, the difference between 
Figures 10a and 10b is small. For the compensated void, there is no perceptible outbound 
shock, and the original density "bump" in the void wall is stretched and damped out. In 
fact, the fluid initially in the compensated void's wall is moving inward at t ~ 2.1. 
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Up to this point, we have shown the evolution of general relativistic uncompensated 



and compensated voids for the initial velocity profile U = t/cRAv = \J2GnM/R (or 
r{ti,R) = 1), where the initial velocity per shell just balances gravity. If the wall ini- 
tially has an outward peculiar velocity (as expected from first-order inflation, for exam- 
ple), then f/(ti,-Rwaii(ii)) > f/GRAv(ti,-Rwaii(^i)) (or r(ti,i?„aii(ti)) > 1), as discussed in 
II-D. In Figure 11, we show inp, M, T and $ versus R R{t\, Rb) / R{t, Rb) for a com- 



pensated superhorizon-sized void with U/c = RJStiGnp/^ (Eqn(|3.5|)), RwaiiiU) = 500, 



Ai?waii(ti) = 10., fc2 = 4, eituRs) = 10^ AR{ti) = 2.5, / = .005 and a = 10"^°. 
The voids are superhorizon-sized: R^^\\{ti) / H {ti, Rb)~^ = 250. In addition, Rjgit) = 
601.2, 783.9, and 911.8, and we show the configurations at ti and at times t = 1.7 and 
t = 2.3. The void wall is initially moving outward with a very large peculiar velocity, 
since r(ti,i?„,u(ti)) ^ 900. 

We find the very interesting result that even though the wall has a large outward 
peculiar velocity, the inner part of the void still collapses. This is because the fluid near 
the base of the void wall gets accelerated into the void right away, pulling adjacent fluid 
with it. It is true however, that at t ~ 2 the density profile looks different than that 
from Figure 10b. This is because the fluid in the void wall takes more time to lose its 
outward velocity. Since the numerical collapse time is Ate = 1-3 from Figure 11, we find 
that the extra time taken for the void to collapse is approximately 1.3 — 1.1 ^ .2, which 
is still less than the outside Hubble time. (This follows qualitatively from our discussion 
in Section II-D, where it was argued that F would decrease very quickly at the inner 
edge of the wall, since F oc — F^p'). These new results show that the collapse time of a 
superhorizon-sized void with a large outward peculiar wall velocity can be of order the 
F^*-crossing time. Since the minimum thermalization and homogenization time is the 
I'^^-crossing time, the time for thermalization and homogenization of this void may be 
short. 

In Figure 12, we show the void radius versus cosmic time for a void with initial 
size -Rwaii(^i) = cW^^H~^^{ti) . If the temperature outside the void initially is Tout(ti) = 
lO^^GeV, then the initial time is t[ = 10~^^ seconds.^! And because recombination oc- 
curs at tree — 10^^ seconds, trcc/^i = 10^^, which is near where the dashed lines in- 
tersect at the top of Figure 12. Therefore, the plot consists the radiation dominated 



epoch in the early universe, where the scale factor outside the void is a{t) = a{t{)Jt/t\ 



Because R{ti, Rb) / Rit, Rb) = \/ti/t, any comoving point (defined as oc Jt/ti, not 



r = const) is a vertical line in this plot. The dashed line labeled rcM shows the 
void size if it were comoving with spacetime, as previously suggested (see IV). The 
Hubble radius outside the void is ciJ^t = -Rhor = 2ct, and is shown as the dashed 
line labeled by thor- As can be seen, thor and vcm intersect at t/ti ~ 10^^ (or 
10~^^ seconds), shortly after recombination. (This also follows from Eqn (|4.2|) , since 
Ate ^ i/o"ut(^i)(c"^i?waii(ti)/i^o"ut(^i))^ - 10^^). If a superhorizon-sized void were to con- 
formally expand with spacetime, then this void would just barely be around to distort 
the microwave background at recombination.! 
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The dash-dot hnes show the position of the I'^^-crossing photons (i.e. the inner void 
wall) from Eqn (|48|) . We show $in(ti) = 5 x 10^\ 10^^ and 10^° (since we require $in(ii) ^ 
VlO^^ (from Eqn (|4.12| ))), with the predicted I'^^-crossing times of tc/^i — AtcAi = 4x10^^, 
10^^ and 10^ respectively. If At^/ H^^^^{t{) > 1, the radius of the inner edge of the wall 
is constant {R ~ const) until time t t^. — ti{^i^{ti)R^^\\{ti) / H~^^{t\))'^ , at which 
point R rapidly decreases to zero (see Eqn (|4.8|) ). If the reheat temperature is Trh = 
Tontiti) = lO^^GeV, then the void with $in(ti) = Tout (ti) /Tin (ti) = 10^° may thermalize 
by t ^ 10^^^ seconds (or at T ~ T^nJti/t ~ lO^^GeV), far before recombination or 



-RH1 

nucleosynthesis. (Note that this value of $in(ii) corresponds to an initial void temperature 
of T^ti) = TRH$in(ti)-' = lO-^GeV = IkeV). 

In conclusion, we have seen that a compensated or uncompensated superhorizon-sized 
void collapses as fluid in the wall diffuses into the void at the speed of light. Thus, the 
collapse time is found to be approximately equal to the I'^^-crossing time calculated in 
Section IV, which is much smaller than previously thought. 

VIII. Discussion 

In this paper, the evolution of a superhorizon-sized void embedded in a Friedmann- 
Robertson- Walker (FRW) universe was studied by numerically evolving the spherically 
symmetric general relativistic equations in the Lagrangian gauge and synchronous co- 
ordinates. (A superhorizon-sized void is a void larger than the Hubble radius outside 
the void: c~^i?waii(^i)/ H~^^{ti) > 1, where ti is the initial cosmic time, cif^t(ti) is the 
Hubble radius outside the void and -Rwaii(^i) is the radius of the void). The particles are 
assumed to be in local thermal equilibrium so that they can be described as a fluid, and 
the perfect fluid equation of state is chosen. The inside of the void is initially chosen 
to be homogeneous and non-expanding. We are particularly interested in the evolution 
of voids with T//i 3> 1, where T is the local temperature and is the mass of a fluid 
particle. This corresponds to a radiation-dominated period in the early universe. We 
find that for p = p/3, general relativistic voids collapse via a shock propagating inward 
from the void wall. The shock is formed from the steep pressure gradient at the inner 
edge of the wall. It moves at approximately the speed of light; this is not surprising, 
since the speed of sound is .57c. Thus the void collapse time (the time taken for the 
shock to reach the origin) roughly equals the photon f^^-crossing time (the time taken 
for a photon initially at the inner wall edge to reach the origin). At the time this shock 
reaches the origin, much of the fluid in the original wall area is moving toward the origin 
behind the shock. The energy density of this fluid (i.e. the fluid behind the shock) is 
less than that outside the void. At the collapse time then, only part of the void has been 
filled in with fluid. In particular, as the initial energy density inside the void decreases, 
the fraction of the void filled in at the collapse time decreases. 

Because the shock moves inward at roughly the speed of light, we can calculate the 
approximate collapse time. For p = p/3, the P*-crossing time is 
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for i?„all(ti)/i^out(^i) ^ \/pout(ti)/Pui(ti), where pin(ti) and pout(ii) are the initial fluid en- 
ergy densities inside and outside the void, respectively, and p = T'^. If c^^i?waii(^i)/-f^cmt(^i) < 
(Pout(^i)/Pin(^i))^'^'^, then /S.tc/ H~^^_{t]) < 1; the f^^-crossing time is less than or compara- 
ble to the initial Hubble time outside the void. In fact, as the density inside the void 
approaches zero, the collapse time goes to zero! 

It may seem contradictory that the I'^^-crossing time for a superhorizon- sized lelaiivis- 
tic void can be less than the outside Hubble time (i.e. or that light travels comparatively 
farther inside a void than outside). There are several points to make about this. First, if 
we choose Eulerian instead of Lagrangian synchronous coordinates to evolve special rela- 
tivistic voids {Gn = 0), then the I'^^-crossing time is not fast. This is because spacetime 
is sliced differently in time inside the void. Since a similar I'^^-crossing time is obtained 
for general relativistic voids, it might be argued that the quick collapse time is a fig- 
ment of the Lagrangian gauge used. However because the FRW metric has a coordinate 
singularity at the Hubble radius when transformed to Eulerian synchronous coordinates, 
superhorizon-sized voids cannot be evolved in these coordinates. Any prior intuition from 
special relativistic voids in these coordinates therefore, must be carefully applied. 

Second, calling a void "superhorizon-sized" is deceptive. Although the f^^-crossing 
time for a superhorizon-sized perturbation in the FRW universe is greater than the outside 
Hubble time, this does not imply that the same is true for a superhorizon-sized void. This 
is because the void is not a small perturbation in general (i.e. (pout — Pin)/Pout — 1). If 
the void's size is defined to be the spacelike circumferential radius relative to the outside 
Hubble radius, then the void is superhorizon-sized if -Rwaii > cH~^^. However, if the 
voids's size is defined to be its radius relative to the Hubble radius inside the void, then 
its size is smaller since c~^R^sii(ti)/ H-~^(ti) = c~^R^aii{ti)/ H~^t{U)\/ Pm(ti)/ Poutiti)] it is 

subhorizon- sized if c~^i?„aii(^i)/-f^<^t(^i) < \/pout(ii)/Pin(ii) ' (We call this latter size the 
"radial size" , because it is the time taken for a photon to cross a static void multiplied 
by c). This is because the Hubble radius is much larger inside than outside the void. 
Since the void is an underdense region, we expect the opposite relative size problem to 
occur for overdense regions where the gravitational potential is large, not small. As an 
example, a black hole's circumferential size is small (or at least finite), while its "radial 
size" is infinite. 

The original motivation for studying this problem was to determine the evolution of 
voids formed from first-order (e.g. extended) inflation. These voids are compensated in 
energy density, and the walls initially have large outward peculiar velocities.tJ Previous 
authors estimated the minimum homogenization time by determining the l^*-crossing 
time. They found it to bJii 

This estimate was based on the assumption that a superhorizon-sized void would con- 
formally expand with spacetime. The present work throws considerable doubt on this 
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assumption and its implications (i.e. a long thermalization time for superhorizon-sized 
voids) for several reasons. First, the l^*-crossing time is much shorter than this estimate. 
For Ate/ H~^^_(t{) > I5 this estimate is off by the factor \J Pmiti) / Pontiti) (see Eqn( p.l| )). 
Since the minimum thermalization time is the I'^^-crossing time, this implies that thermal- 
ization and homogenization may happen much quicker than suggested. This would have 
profound effects upon our understanding of the evolution of superhorizon-sized voids in 
the early universe. Second, the inner edge of the void collapses (not expands) in roughly 
the l***-crossing time, even (apparently) for voids with large outward peculiar veloci- 
ties. This collapse occurs because the wall decelerates from the positive wall pressure 
{U (X —p' < 0). If the wall pressure were negative instead (e.g. the wall pressure of a 
vacuum bubble), the wall would accelerate outward. Thus, fluid from at least part of the 
wall rushes into the void, partially "filling it" at the I'^^-crossing time. Because the light 
travel distance is still large there, the void may still fill up in a relatively short amount 
of time. Third, the void wall can cease expanding by many mechanisms. Because fluid 
rushes into the void, large velocity gradients are formed in the original wall area which 
"pull" on the wall and can slow it down. In addition, when the peculiar wall velocity is 
much larger than the gravitational velocity, the deceleration of the wall is proportional to 
the velocity squared which can be enormous. Finally, an outgoing wave will be damped 
and slowed down in any case, by virtue of the volume effect caused by the spherical 
geometry. Thus even if the wall initially expands outward conformally (or faster) , it may 
stop doing so rather quickly, as the simulations appear to suggest. In any case, even if 
the outer part of the wall can carry out some mass for a long period of time, the interior 
region will continue to fill in and homogenize. Since overdensities will most likely form 
near the origin after collapse, perturbations, small waves and/or overdensities will most 
likely be around at recombination. These perturbations could have interesting ampli- 
tudes, and therefore could alter the standard Harrison-Zeldovich density spectrum. In 
any case, because the qualitative void evolution picture previously suggested is at least 
partially incorrect, it is reasonable to expect the estimated thermalization time to be 
erroneous; the "big-bubble problem" might not be a problem after all. 

It is important to point out that in order not to distort the microwave background, 
the thermalization time only needs to decrease somewhat from the estimate given by 
Eqn (|8.2|) . The largest superhorizon-sized void from minimalist inflation is roughly c^^ 
Rwaii{ti) / H~^^-{ti) ~ 10^''. la If the initial time after reheat is ti = 10~^^ seconds (Trh = 
lO^^GeV), then at recombination (tj-ec = 10^^ seconds), this void would not have ther- 
malized according to Eqn( |8.2| ). Using these initial conditions with Eqn( |8.2| ), it is easy 
to see that voids for which c~^-Rwaii(^i)/-f^(^t(^i) ~ 10^^ were traditionally troublesome. 
However, since 10^^ is such an enormous number, a very slight change in the functional 
form for the thermalization time, Atn, can cause thermalization to occur before recom- 
bination. For instance, if we suppose that the thermalization time can be written as 
Atn ~ c~^H~^^{ti){R^i^ii{ti)/ H~^^(ti)yp, then an acceptable range for p is p < 1.6, which is 
not dramatically different from 2. The point is that the "big-bubble problem" can only 
be resolved by accurate knowledge of the thermalization time, because estimates off by 
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a somewhat small amount can lead to erroneous conclusions. 

Of the many questions which remain unanswered, the most important is to deter- 
mine what happens to a superhorizon-sized void after it collapses. If it thermalizes and 
homogenizes, at what time does this happen? As mentioned previously, at the time of 
l*^*-crossing, the energy density within the former void region is greater than Pin(ti), but 
is still less than Pout(^)- In addition, the bulk motion of much of the original wall fluid 
is inward so that spacetime is collapsing and not expanding there; thermalization and 
homogenization has not occurred. However, since the "void" is still relatively underdense 
at this time, the light travel distance will be comparatively large, thereby allowing for 
the possibility of relatively fast thermalization and homogenization. We have found that 
special relativistic voids eventually nearly homogenize after creating an overdense region 
at the origin. In this case, all but the fluid near the origin becomes roughly homogeneous 
fairly quickly. The fluid near the origin however, is overdense for a relatively long time. 
If general relativistic voids behave similarly, then after collapsing, an overdense region 
would form at the origin. This would take a relatively long amount of time to diffuse 
away, and would eventually result in a nearly homogeneous and isotropic universe. In any 
case, since this overdense region would form very quickly, it is unlikely that empty voids 
would be around any time near recombination, as previously suggestedSii Instead, 
overdensities, perturbations and small waves with strange fluid velocities may be. Future 
work will study the consequent evolution of a general relativistic void after collapse,!! 
as well as the evolution of a superhorizon-sized void from first-order (e.g. extended) 
infiationS 
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Appendix A: Transformation of the FRW metric to Eulerian Coordinates 

As discussed in Section IIA., the metric in Eq. ( p.lD can be subjected to transforma- 
tions of the type t = fi{t',r') and r = f2{t',r') without altering its spherical symmetry. 
If we set = k{t,r), then we can rewrite Eq. (|2.1|) as 

ds^ = -{N^ ^ S^)dt^ + 2Sdtdr + h^dr^ + r^dQ'^, (A.l) 

where we have introduced the lapse and shift functions as N(t, r) and S(t, r), respectively. 
This is the Eulerian gauge with synchronous coordinates if we choose S = 0. As we shall 
see, when the FRW spatially flat metric is transformed in this manner, the resulting 
metric has a coordinate singularity at the horizon. Thus, the study of superhorizon-sized 
voids is not possible in these coordinates. 

Consider the FRW k = metric given by Eqn ( p.3|) . Letting f = ra{t), we can rewrite 



29 



this as 

ds^ = - (^1 - f2 (^^) dt^ - 2r-Jirdt + dr^ + rHn\ (A.2) 

where we set c = 1 everywhere in this appendix. This is the Eulerian asynchronous 
form of the FRW metric. We notice that this takes the form of Eq. with = 1, 

S = rd/a and h = 1. Since a{t) oc t^, where ^ = 1/2 and 2/3 in radiation (p = p/3) 
and matter {p = 0) dominated universes respectively, we can ehminate the drdt term by 
defining 

t = ^r/g{i,r) (A.3) 

and ensuring that the condition 

^ + {1-0 9^ 
dr ^f(l — g"^ 

hold. Defining k = (1 — 0/^ a = 1/[2(1 — ^)], this last expression can be integrated 
to obtain 

h(ty=g/il + Kg'r, (A.5) 
where h{t) is an arbitrary (integration) function of t only. The metric then becomes 



9 'Z^r'i (A.4) 




ds^ = _l'^\ ^'^^y ) dP + ^-dr^ + r^dn\ (A.6) 



' l +J^g^f^,r2 

1-9' ' 1-9'' 

This metric has a coordinate singularity at (7 = 1, which from Eqn( |A.3|) occurs at 
r = t/a(t) But this is the comoving radius of the Hubble radius (rnoR), since 

H^^ = t/^ = rHORO'(^)- Therefore, we have shown that for the FRW models expressed in 
Eulerian, synchronous coordinates, a coordinate singularity occurs at the Hubble radius. 

As an interesting example we take ^ = 1/2, which corresponds to radiation-domination 
{p = P/3)- In addition, we choose h = ^/i. Using Eqs. ( |y4.5| ) and ( |A.3| ) together with 
r = r a{t) = r a{t\){t/ti)^/'^ , we find that a = K = l, t = t + {a{ti)ry / (Ati) and 

g = {1 — \Jl — A[r /ty) / {r /t) . In addition, the metric simplifies to ds^ = — g'^)~^d'P + 
(1 — g'^)~^dr^ + r^dfl"^. Note that even though the energy density is spatially constant 
on a hypersurface t = constant, on a hypersurface of constant t, p cc t~^ = 4g'^r~'^ = 



4p(l — y 1 — A{r/ty )^/f^ which is not spatially constant. 

Appendix B: Differencing of the Equations 

In this appendix, we list the expressions as they are differenced in the code. We use 
the shorthand notation AG = — Gj or AG = Gj — Gl_i for forward or backward 
differencing, respectively. The equations involving spatial derivatives are 

n^A{R'^U) 

(B.l) 



n 



R^AR 

M{p + Q)A{R^U) 

Tr'^Ar 
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We note that there are many other schemes which could potentially be used. (For ex- 
ample, one could write instead of r^Ar in the expression for e ). Although these 
variations can produce and propagate the shock and rarefaction waves, they can not han- 
dle the bounce of the shock at the origin. The only set of difference equations which we 
have found to do this are given above. For the Tolman-Bondi pressureless dust models 
(Section V), it is found that more accurate solutions are obtained near the origin when 
differencing the n equation as follows: n = — (n^$A(i?^f/))/(rr^Ar). The Tolman-Bondi 
figures shown in this paper, however, use the difference equations as written in Eqns( [BTT|) . 
Appendix C: General Relativistic Jump Conditions 

In this appendix, we derive the jump conditions for general relativistic shocks. This 
closely follows the derivation by May and White (1965)lij. These conditions will then be 
used to derive a more general artificial viscosity expression. In addition, the shock jump 
conditions for ultra special relativistic shocks are derived. 

Let the variables a and b represent labels for comoving observers ahead and behind 
the shock, respectively. If each observer measures the invariant interval separating the 
same two events on the world surface of a shock, using Eqn(|2.2|), we obtain 



[ds'] = [-c'^'dr + A'dr' + R'dn'] = 0, (C.l) 

where [G] = Ga — Gb- Because the shocks are radial, we can choose the two events to 
have dr = dt = but dQ ^ 0. Therefore, R is continuous across the shock: [R] = 0. 
Now consider two events with dr ^ and dt ^ 0. Then from Eqn (|C.l|) , 

[c^$2 - M^/n^] = 0, (C.2) 

where is the position of the shock in comoving coordinates, drg/ dt is the shock "speed" , 
and Ms = f {drJdt)/{A'KR'^) from Eqns(|l|) and (^^. This is the first of the jump 
condition equations. 

It can be shown that thejump conditions for the Schwarzschild metric (but not for 
the comoving metric) are Ha 

{T/dg/dx>'] = 0, (C.3) 

where g is the equation for the world surface of the shock. Therefore we must first find 
the jump conditions in the Schwarzschild frame and then transform back to the comoving 
frame. The Schwarzschild metric is ds"^ = —c^A^dT"^ + B^dR^ + R^dfl^, where R is now 
the "Eulerian" coordinate radius. Since i? is a coordinate, [R] = 0, and using the same 
argument as above for the two observers, 

[c^A^ - B^S^] = 0, (C.4) 

where S is the shock "speed" in this frame: cS = dRs/dT. For the Schwarzschild metric, 
the solution for the metrics functions are well known: B^ = 1 — 2MG/{Rc^) = A~'^ 
where M{R, T) = Attc-'^ pR^dR. As long as p is not infinite in the shock, the mass is 
continuous across the shock, [M] = 0. Then [B] = [A] = 0. And using Eqn ( |U.4| ), we 
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find that [S] = 0. We conclude that in the Schwarzschild frame, the metric components 
and the shock "speed" are continuous across the shock: [A] = [B] = [R] = [S] = 0. 

We now derive the jump conditions in the Schwarzschild frame. If a shock is located at 
position Rs(T) at time T, the equation for the world surface of a shock is g = Rs{T) —R = 
0. In addition, the perfect fluid stress-energy tensor in the Schwarzschild frame^ is 



—2/ I \ I la t\ , I ip\ I la rX , / luX 

c {p + p)g^xt^^u +pg^^g =nwg^^u^u +pg^),g , 



(C.5) 



and the conservation of mass equation is [nu'^dg /dx^] = 0. 00 Using Eqns ( |C.3| ), the 
junction conditions become 
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{-A^u'^fnw + p)S + nwA\'^u"' 
SnwE^'^u"" - {nwB\u'y + p) 
nSu'^ - nu'^ 





0. 



(C.6) 
(C.7) 
(C.8) 



We would like to express these equations in terms of the comoving metric functions. 
First, we can rewrite the shock velocity in the comoving frame as cS = {Rt + RtTs) / {Tt + 
Trfs). In addition, we can relate the metric functions from the comoving frame to those in 
the Schwarzschild frame by g'^" = {dx'^ /dx'^) {dx'" /dx^) g"^. We then obtain (cA)"^ = 
(c$)-2t2-A-2t2, = -{c^y^Rl + k'^Rl and {c^Y^TtRt = A-^T.i?,. In addition, 
the "mass" M(r, t) contained within comoving radius r is continuous across the shock, 
since [M] = [nR'^{l + e/c^)]dR ^^^q 0- Using Eqn ( gl2D , we find that [T^ - U^/c^] = 0. 

We can use these relations to calculate the 4- velocity as measured in the Schwarzschild 
frame. Since the 4-velocity in the comoving frame is = (— c$-^, 0, 0, 0), and the ve- 
locity transforms as u'^ = {dx'^/dx") u" , u'^ = -^-\cf , R, 0,0), where f = dT/dt 
and R = dR/dt. Since U = -R/$, and using the fact that the 4-velocity is normal- 
ized to u'^'^u'x 



— c 



c A {u ) — B [u ) , we can rewrite the 4-velocity as u 



iX 



-(A-V52[/2 + c2,f/,0,0). In the special relativistic limit (G^v 
and this becomes 



u 



iX 



-{{U' + c'f/\U,0,0) = -{cT,Tv,0,0), 



0), A = 1 and B = 1 



(C.9) 



where we have defined v = U /T so that F = (y 1 — (f/c)2)~^. In the special relativistic 
limit, then, v is the fluid particle's radial velocity, F is usual gamma-factor (i.e. energy 
per particle mass), and U is the particle momentum per particle mass. 

We can now rewrite Eqn ( |(J.6|) - ([(T8D and [S] = in terms of the comoving metric 
functions using the fact that T^/T^ = c~2[/A/(F$) and [r^] = 0. We obtain 



[UMs/{n^) + <I>F] 
[c2M,F(1 + e/c^) - pU^ 
[MsU{l + t/^)-pT^ 






0. 



(C.IO) 
(C.ll) 
(C.12) 
(C.13) 



°We denote quantities in the Schwarzschild frame by primes. 
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Eqn ( |C.2| ) and Eqns ( |C.10| )-( |C.13| ) make up the required shock conditions. One of them 
however, is redundant. It can be shown that in the nonrelativistic hmit, the above 
conditions reduce to the Lagrangian shock jump conditions given in Ref 

The case considered by May and White (1967) is the penetration of a shock into a 
non-relativistic medium. We set Gn = and ea = Pa = Ua = ^ and take e = p/[{pf — l)n]. 
Using Eqns ( |C.11|) and ( |C.12|) with the fact that U'^ = — 1, the energy behind the 



shock is found to be = Ff, — 1. Using Eqns (|C.12|) and ( |C.13| ), we ehminate and 
find that t] = rib/na = {1 + Ffe7)/(7 — 1). From Eqn ( |U.12|) , we find that 

M = (7 - lKeb$fe/f/b, (C.14) 

and plugging this into Eqn (|C.2| ), we determine that $b = 1/(1 + 7(F;, — 1)). Finally, we 



plug this into Eqn (|Cl4l ) and find M, = cn„ J{Tb - l)/(Fb + 1) (l + Ffc7)/[l+7(Fb-l)] 



(It turns out that Eqn ( |C.10| ) gives no new information). Thus, all quantities behind the 
shock can be expressed in terms of Ff,, the "strength" of the shock. 

How does this affect the artificial viscosity needed to smooth the shock? Von- 
Neumann's viscosityll is Q = k'^niJJ'Ydr'^ for f/' < and otherwise, with the shock 
being spread over roughly k"^ grid points. Over the k"^ points, Q ~ n{Ub — UaY = nll^ = 
c^n{Tl - 1) = c^niTh - l)(Ffc + 1). Using the fact that e^/c^ = Ff, - 1, Q ~ nehTb ~ pT 
which is greater than the pressure p by the factor of F > 1. Since Q should be of order the 
pressure in the shock region, we divide by F: Q = k'^n{U'Ydr'^ /T for [/' < and Q = 
otherwise. In the nonrelativistic limit, F = 1 and Von-Neumann's artificial viscosity is 
obtained. 

In this paper, we are more interested in the case where the fluid on both sides of 
the shock is ultrarelativistic. We first set Gat = and Ua = 0. Then Eqn( |C.2| ) and 
Eqns ( P^^) - (|C . 1 3| ) can be manipulated to yield 



Ub = JTI-I (C.15) 



M = aiaTi^Tl-l^J {VbV~l) (C.16) 

r/ = Tb{l + iebl^)-{l + eJ^) / [eblc\^-l)] (C.17) 

eb/c^ = -l + Vb{l + ta/c^) + {i-l){Vb-ri-^)ta/c^ (C.18) 

= \T^~Tb\^a / {Tbr]-l). (C.19) 

(It can be shown that these equations reduce correctly in the — > limit). Now, if we 
set 7 = 4/3 and take the ultrarelativistic limit, e^/c^ > ea/c^ » 1, Eqn ( p.l7| ) becomes 



rj = 4Fft{l — 3ea/{4Tbeb)} — 4Tb. Then the equations can be approximately solved to give 

eb ~ 4eaTb{l - (4F,)-2) / 3, (C.20) 

which to lowest order is eb — 4eaFfe/3. Therefore, eb depends not only on F;, but also on 
ea- We generalize the artificial viscosity to 

Q = en (l + e/ (Fc^)) {U'fdr^/ F for U' < 

g = otherwise. (C.21) 
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For a strong shock then, Q ~ neU"^ / (F^c^) ne p, as desired. And when —>■ oo, this 
becomes VonNeumann's original expression. Note that because this was derived in the 
hmit of strong shocks, it does not work as well for weaker special or general relativistic 
shocks. ii 

We now write down the exact solution of Eqns ( p.l5| )-( [C.19| ) for conditions behind the 



shock in terms of Tf, only and quantities in front of the shock. Because p = ne in this limit, 
we can rewrite Eqns ( P.17|) and (|C.l^ ) as 77 = 4:pbTb/{'iPa + Pb) and r] = {Spb+ pa) / {4:paTh) 
respectively. Setting them equal, we can solve for p^. Then, all others quantities are 
determined. They are 



Ub = (C.22) 
Pb - 10 



Pa 6 



1 + a/1 -36 (IGP-IO)- 



(C.23) 



4rbY'3pb/pa 
nh 4rf, 



M = + ^-'^"(3p,/p. + 1) ^^24) 



(C.25) 
(C.26) 



Ha 3(pb/pa)~^+l 
^ + Pb/Pa 



en 4r 



b 



= ^■t„^ =' + y°>" fl... (C.27) 

efe 4rb 

We note from Eqn ( |(I2^ ) that [^w] = 0. 

Appendix D: Nonrelativistic Shock Tube Problem 

In this section, we sketch the derivation of the slab shock tube solution for non- 
relativistic fluids.E3 We start with a fluid in which the pressure p and specific volume 
V = 1/n (where n is the number density) to the left and right of Xq are pi, V\ and p^-, V2 
resectively — thus, p and V are discontinuous across x = xq. We assume here that our 
coordinates are oriented so that p\ < p2 and Vi > V2. In addition, the fluid is initially at 
rest: f 1 = f 2 = 0. Because there is no scale to the problem, the solution is a function of 
x/t only, where x is the Eulerian coordinate. At time t later, there are 4 regions. Region 
1 and 2 are at rest with p = pi,V = Vi and v = Vi, and ^ = ^25^ = ^2 and v = V2 
respectively. Separating region 1 and 3 is a shock wave, with position Xd- Region 3 
contains the fluid behind the shock, and region 3' is wedged between regions 3 and 2 and 
contains the rarefaction wave. The boundary between region 3 and 3' is called a contact 
discontinuity and is at position Xc- The velocity and pressure are constant across this 
boundary: py = p^ and = ^3. The number density, however, is not. Finally the 
boundary between region 3' and region 2 is located at Xa- We will calculate p, V and v 
for regions 3 and 3' at time t, assuming that e = pV/{'j — 1). We will solve for the 
following unknowns: ps, V3, n^, ny and p, ra, and v in the rarefaction wave. 

Across a shock front, t;i --Us = J {p^ - VijiVi - ^3) and ei-e3+|(Vi-V3)(pi+p3) = 0. 
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Using e = pV/ (7 — 1), they become 



^3 = Vi[{l + + (7 - l)P3]/[(7 - + (7 + 1)P3] (D.l) 
= -{P3 - Pi)pVi/[{^ - + (7 + 1)P3] (D.2) 

The similarity solution for a rarefaction wave is 

p = p2[l-(7-l)l^l/(2c2)]^/(^-^) (D.3) 

^ = V2{p2/p)y/^ (D.4) 

= 2(c-C2)/(7-l) (D.5) 

|.;| = 2(c2-x/t)/(7 + l), (D.6) 

where the local speed of sound is given by Cj = y/'jpiVi. We can equate the speed of the 
fluid V3 to the fluid velocity in the rarefaction wave between regions 3 and 3'. We find 

v,y = 2C2[(P2M)^'-"^/^'"^ - l]/(7 - 1) (D.7) 
Combining Eqns (p.2| ) and (p.7|) , we find that 

P3=Pi + J^^V(7 - l)Pi + (7 + 1)P3 [1 - (P3M)(^-^)('^T . (D.8) 



This equation can now be solved numerically for ^3. We can then determine ^3, n^, and ny 
using Eqns (p.2| ), ( p.l|) and (|D.5| ). In addition, we find p, n and v from Eqns (p.3| ), 
( p.4| ) and ( |D.5| ). in the rarefaction wave. 

The final step to the solution is determining the location of the boundaries separating 
these regions. Since the rarefaction wave is moving at the sound speed in region 2, 
Xa = tc2. Since the velocity at 6 is f = v^' from Eqn ( p.7| ), we then use Eqn ( |D.6|) 
to find its position: Xb = tc2[l — (7 + l)/(7 — 1) [1 — iP2/P3)^^~"'^^'^'^^]]- The contact 
discontinuity moves with the fluid and its position is thus Xc = V3 t. Finally, the shock 
position can be determined by transforming to a frame in which the shock is constant 
and using mass conservation. We then transform back to find the the velocity of the 
shock to be Vg = (ns/ (723 — ni))v3 so that x^ = Vgt. 
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Fig. 1: The evolution of a compensated, pressureless void with radius -R„aii(ti) = 



333c/foul(ti) and for velocity U{ti,R) = J2GnM/R. The energy density function of 




FIGURE CAPTIONS 
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RR{ti, RB)/R{t, Rb) — R{ti/ty/^ is shown at the initial time U — 1 and at times t — 10, 
100 and 300. 

Fig. 2: Same as for Figure 1 but with velocity U{ti, R) = c~^R^JStiGnp/'^- A dense, 
thin, outward-moving shell is formed. Shell-crossing occurs t — 1.048. 

Fig. 3a: The evolution of a pressureless, uncompensated void with radius -Rwaii(^i) — 
.0067cifout(^i) s-iid velocity [/(tj, R) — c'^R^JSirGisfp/S. Artificial viscosity is used to pre- 
vent shell-crossing. The energy density as a function of RR{ti, RB)/R{t, Rb) — R{ti/ty^^ 
is shown at the initial time U — l and at times t — 50, 100, 400 and 1000. 

Fig. 3b: Same as Figure 3a but for a compensated void. 

Fig. 3c: The radius of the shell versus time for the results pictured in Figures 3a and 
3b. The predicted self-similar solutions arc overlaid as a dotted and dashed line for an 
uncompensated (i?sheii oc t^^^) and compensated (-Rgheii c< t^^^) void respectively. 

Fig. 4: The relative errror in the energy density versus RR(ti, Rb)/ R{t, Rb) — 

R^Jti/t for an initially homogeneous and isotropic FRW universe with p = p/3 at time 

At/ti = .1. The results are for / = .01, .005, and .0025 shown with solid, dotted and 
dashed lines, respectively. 

Fig. 5: The nonrelativistic shock tube for -Rwaii(ii) = 20, Ai?„aii(^i) = -01 and 
Gjv = 0. The dashed lines represent the exact slab similarity solution. 

Fig. 6a: The nonrelativistic evolution of an uncompensated void with Gn = 0. The 
dashed line is the initial time ti = 1, and the triangles with connecting lines is at time 
t = 2.5 after the inbound shock has bounced and the fluid has begun to settle down. 

Fig. 6b: Same as Figure 6a but for a compensated void. 

Fig. 7a: Collapse of a special relativistic, uncompensated void for i?waii(^i) = 1; 
Ai?waii(ii) = -02, T/fjL > 1, Gat = 0, U{ti,R) = and a = 10"^ Shown is M, F, Anp 
and $ versus R for the initial time ti — 1 and for times t — 1.04 and 1.065 (before the 
shock reaches the origin). 

Fig. 7b: We plot Airp versus R for special relativistic voids with initial varying wall 
thicknesses at the initial time (unlabeled) and at later times t = 1.04 and 1.07, 1.08 or 
1.065. The parameters are the same as in Figure 7a but with Ai?waii(^i) — -02, .04, .06 
and .1. 

Fig. 7c: Collapse, thermalization and homogenization of a special relativistic void. 
The parameters are the same as in Figure 7a. Shown is the initial time t^ and times 
t = 1.04 (void collapsing) and t = 3.0 (after collapse). 



38 



Fig. 8: Collapse of a general relativistic, uncompensated void for -Rwaii(^i) = 1, 
A-Rwaii(^i) = -02, T/n > 1, U{ti,R) = yj2GNM/R and a = lO"'' at the initial time 
ti — 1 and at times t — 1.04 and 1.065. 

Figs. 9a,b,c: Collapse of general relativistic, uncompensated voids for -Rwaii(^i)/ 
H-^^{U) = 25, Ai?^aii(ii) = 1, T/n > 1, U{ti,R) = ^2GnM/R. Figures 9a, 9b and 9c 

show the pressure versus RR{ti, RB)/R{t, Rb) — R\JuJt for a — 10~^, a — 10~^ and 
a — 10"-^°, respectively. 

Figs. 10a,b: Collapse of general relativistic, uncompensated and compensated void 
in 10a and b , respectively. i?waii(^i)/^o"ut(^i) = 250, Ai?^aii(^i) = 10, T/n > 1, U{t;,R) = 
^2GnM/R and a — 10~^° for both. Shown is the pressure versus RR{ti, RB)/R{t, Rb) — 

R\fkft. 

Fig. 11: Collapse of a general relativistic, compensated void for i?„aii(^i)/-f^out(^i) — 
250, Ai?^aii(ti) = 10., T//X > 1, U{U,R) = c-^R^StiGnp/S and a = lQ-^°. Shown is 
the pressure versus RR{ti, Rb)/ R{t, Rb) — R^U/t. 

Fig. 12: l^*-crossing times versus R^U/t for general relativistic voids with initial 

radius -Rwaii(^i) = cl0^^i7Jut(^i) outside temperature Tout = 2rh, where Trh is the 
reheat temperature. The dash-dot lines are for $in(ti) = 7'out(*i)/^in(^i) = 5 x 10^^, 10^^ 
and 10^° and the temperature at which each void potentially thermalizes is T = 500, 
10^ and lO^^GeV respectively. We also plot (long dashes) the evolution of a conformally 
stretched perturbation with the same initial size, rem, and the Hubble radius (short 
dashes) outside the void, thor- H ^rh = lO^^GeV, then recombination occurs at t/ti = 
10^^ (or at temperature Tree — 3 x 10~^GeV), where rem and thor intersect. 
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